Numerical Solution of Ordinary 
Differential Equations 

E. Siili 
August 12, 2010 



1 



Contents 



1 Picard's theorem 1 

2 One-step methods 4 

2.1 Euler's method and its relatives: the ^-method 4 

2.2 Error analysis of the 0-method 7 

2.3 General explicit one-step method 9 

2.4 Runge-Kutta methods 13 

2.5 Absolute stability of Runge-Kutta methods 19 

3 Linear multi-step methods 21 

3.1 Construction of linear multi-step methods 22 

3.2 Zero-stability 24 

3.3 Consistency 26 

3.4 Convergence 29 

3.4.1 Necessary conditions for convergence 30 

3.4.2 Sufficient conditions for convergence 33 

3.5 Maximum order of a zero-stable linear multi-step method 37 

3.6 Absolute stability of linear multistep methods 43 

3.7 General methods for locating the interval of absolute stability 45 

3.7.1 The Schur criterion 45 

3.7.2 The Routh-Hurwitz criterion 46 

3.8 Predictor-corrector methods 48 

3.8.1 Absolute stability of predictor-corrector methods 50 

3.8.2 The accuracy of predictor-corrector methods 52 

4 Stiff problems 53 

4.1 Stability of numerical methods for stiff systems 54 

4.2 Backward differentiation methods for stiff systems 57 

4.3 Gear's method 58 

5 Nonlinear Stability 59 

6 Boundary value problems 62 

6.1 Shooting methods 62 

6.1.1 The method of bisection 63 

6.1.2 The Newton-Raphson method 63 

6.2 Matrix methods 66 

6.2.1 Linear boundary value problem 66 

6.2.2 Nonlinear boundary value problem 69 

6.3 Collocation method 70 



Preface. The purpose of these lecture notes is to provide an introduction to compu- 
tational methods for the approximate solution of ordinary differential equations (ODEs). 
Only minimal prerequisites in differential and integral calculus, differential equation the- 
ory, complex analysis and linear algebra are assumed. The notes focus on the construction 
of numerical algorithms for ODEs and the mathematical analysis of their behaviour, cov- 
ering the material taught in the M.Sc. in Mathematical Modelling and Scientific Compu- 
tation in the eight-lecture course Numerical Solution of Ordinary Differential Equations. 

The notes begin with a study of well-posedness of initial value problems for a first- 
order differential equations and systems of such equations. The basic ideas of discretisation 
and error analysis are then introduced in the case of one-step methods. This is followed 
by an extension of the concepts of stability and accuracy to linear multi-step methods, 
including predictor corrector methods, and a brief excursion into numerical methods for 
stiff systems of ODEs. The final sections are devoted to an overview of classical algorithms 
for the numerical solution of two-point boundary value problems. 

Syllabus. Approximation of initial value problems for ordinary differential equations: 
one-step methods including the explicit and implicit Euler methods, the trapezium rule 
method, and Runge-Kutta methods. Linear multi-step methods: consistency, zero- 
stability and convergence; absolute stability. Predictor-corrector methods. 

Stiffness, stability regions, Gear's methods and their implementation. Nonlinear stability. 

Boundary value problems: shooting methods, matrix methods and collocation. 
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1 Picard's theorem 



Ordinary differential equations frequentfy occur as mathematical models in many branches 
of science, engineering and economy. Unfortunately it is seldom that these equations have 
solutions that can be expressed in closed form, so it is common to seek approximate 
solutions by means of numerical methods; nowadays this can usually be achieved very in- 
expensively to high accuracy and with a reliable bound on the error between the analytical 
solution and its numerical approximation. In this section we shall be concerned with the 
construction and the analysis of numerical methods for first-order differential equations of 
the form 

y' = f{x,y) (1) 

for the real- valued function y of the real variable x, where y' = dy/dx. In order to select 
a particular integral from the infinite family of solution curves that constitute the general 
solution to (1), the differential equation will be considered in tandem with an initial 
condition: given two real numbers Xq and yo, we seek a solution to (1) for x > xq such 
that 

y{xo) = yo ■ (2) 

The differential equation (1) together with the initial condition (2) is called an initial 
value problem. 

In general, even if /(•, •) is a continuous function, there is no guarantee that the 
initial value problem (1-2) possesses a unique solution^. Fortunately, under a further mild 
condition on the function /, the existence and uniqueness of a solution to (1-2) can be 
ensured: the result is encapsulated in the next theorem. 

Theorem 1 (Picard's Theorem^) Suppose that /(•,•) is a continuous function of its 
arguments in a region U of the (x, y) plane which contains the rectangle 

R = {{x, y) : xq<x< Xm, \y - yo\ < Ym} , 

where Xm > xq and Ym > are constants. Suppose also, that there exists a positive 
constant L such that 

\f{x,y)-f{x,z)\<L\y-z\ (3) 
holds whenever {x,y) and {x,z) lie in the rectangle R. Finally, letting 

M = max{\f{x,y)\ : {x,y) G R} , 

suppose that M{Xm — xq) < Ym- Then there exists a unique continuously differentiate 
function x i— > yix), defined on the closed interval [xo,Xm], which satisfies (1) and (2). 

The condition (3) is called a Lipschitz condition^, and L is called the Lipschitz 
constant for /. 

We shall not dwell on the proof of Picard's Theorem; for details, the interested reader 
is referred to any good textbook on the theory of ordinary differential equations, or the 

^Consider, for example, the initial value problem y' = y^^^ , y{0) = 0; this has two solutions: y{x) = 

and y{x) = a;^/27. 

^Emile Picard (1856-1941) 
^Rudolf Lipschitz (1832-1903) 
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lecture notes of P. J. Collins, Differential and Integral Equations, Part I, Mathematical In- 
stitute Oxford, 1988 (reprinted 1990). The essence of the proof is to consider the sequence 
of functions {yn}^=o, defined recursively through what is known as the Picard Iteration: 

yo{x) = yo , 

(4) 

rx 

yn{x) = yo+ n = l,2, 

Jxo 

and show, using the conditions of the theorem, that {yn}^o converges uniformly on the 
interval [xo,Xm] to a function y defined on [xo,Xm] such that 



y{x) = yo+ ^/(e,y(6)d^ 



This then implies that y is continuously differentiable on [ojcXm] and it satisfies the 
differential equation (1) and the initial condition (2). The uniqueness of the solution 
follows from the Lipschitz condition. 

Picard's Theorem has a natural extension to an initial value problem for a system of 
m differential equations of the form 

y' = f(x,y), y(a;o) = yo, (5) 

where yo G R"* and f : [xo,-^jv^] x R"* — R"*. On introducing the Euclidean norm || ■ || 
on R"* by 

(m \ 1/2 

we can state the following result. 

Theorem 2 (Picard's Theorem) Suppose that f(-,-) is a continuous function of its 
arguments in a region U of the {x, y) space R-'^+'" which contains the parallelepiped 

f< = {{x,y) : xo < X < Xm, ||y - yoll < ^m} , 

where Xm > xq and Ym > are constants. Suppose also that there exists a positive 
constant L such that 

||f(ar,y)-f(x,z)|| <L||y-z|| (6) 
holds whenever {x,y) and (x, z) lie in R. Finally, letting 

M = max{||f(x,y)|| : (a;,y) G R} , 

suppose that M{Xm — xq) < Ym- Then there exists a unique continuously differentiable 
function x i— > y{x), defined on the closed interval [xo,Xm], which satisfies (5). 

A sufficient condition for (6) is that f is continuous on R, differentiable at each point 
(x,y) in int(R), the interior of R, and there exists L > such that 



<L for aU (x,y) G int(R) , (7) 
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where di/dy denotes the m x m Jacobi matrix of y G R"* i— f(a;,y) G R"^, and || • || is a 
matrix norm subordinate to the EucUdcan vector norm on R™. Indeed, when (7) holds, 
the Mean Value Theorem implies that (6) is also valid. The converse of this statement is 
not true; for the function f(y) = (|yi|, . . . , \ym\)'^, with xq = and yo = 0, satisfies (6) 
but violates (7) because y t-^ f(y) is not differentiable at the point y = 0. 

As the counter-example in the footnote on page 1 indicates, the expression |y — z| 
in (3) and ||y — z|| in (6) cannot be replaced by expressions of the form \y — z\" and 
||y — z||", respectively, where < a < 1, for otherwise the uniqueness of the solution to 
the corresponding initial value problem cannot be guaranteed. 

We conclude this section by introducing the notion of stability. 

Definition 1 A solution y = v(x) to (5) is said to be stable on the interval [xo,Xm] 
if for every e > there exists S > such that for all z satisfying ||v(xo) — z|| < S the 
solution y = w(a;) to the differential equation y' = f(x,y) satisfying the initial condition 
w(xo) = z is defined for all x G [xqj-'^m] and satisfies ||v(x) — w(x)|| < e for all x in 
[xo,Xm]- 

A solution which is stable on [xo,oo) (i.e. stable on [xo,Xm] for each Xm and with 6 
independent of Xm ) is said to be stable in the sense of Lyapunov. 

Moreover, if 

lim ||v(x) — w(a;)|| = , 

then the solution y = v(a:;) is called asymptotically stable. 

Using this definition, we can state the following theorem. 

Theorem 3 Under the hypotheses of Picard's theorem, the (unique) solution y = v(x) 
to the initial value problem (5) is stable on the interval \xq,Xm\, (where we assume that 
—GO < Xo < Xm < oo). 



Proof: Since 



and 



it follows that 



v{x)=v{xo)+ rmMO)dC 

JXQ 

w(x) = z+ rf(C,w(0)d^, 

Jxn 



||v(x)-w(a;)|| < ||v(xo)-z||+ r||f(e,v(0)-f(^,w(0)||dC 

JxQ 

< ||v(xo)-z||+L ^||v(e)-w(0||d^ (8) 

Jxo 

Now put A{x) = ||v(x) — w(x)|| and a = ||v(xo) — z||; then, (8) can be written as 

A{x) <a + L r A{^) d^ , xo<x<Xm . (9) 

Jxo 

Multiplying (9) by exp(— Lx), we find that 



dx 



Jxo 



< ae--^^ . (10) 
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Integrating the inequality (10), we deduce that 

^-Lx r d^ < y (e--^^° - e"^^) , 
Jxo L ^ ' 

that is 

L r A{0 dC < a (e^(^-^o) - l) . (11) 

Jxo ^ ' 

Now substituting (11) into (9) gives 

^(x) < ae^(^-^°\ xo<x<Xm. (12) 

The imphcation "(9) =^> (12)" is usually referred to as the Gronwall Lemma. Returning 
to our original notation, we conclude from (12) that 

||v(x) - w(x)|| < ||v(xo) - z||e-^(^-^o) , xo<x<Xm . (13) 

Thus, given e > as in Definition 1, we choose 6 = eexp(— L(Xm— a^o)) to deduce stability, 
o 

To conclude this section, we observe that if cither xq = — oo or Xm = +oo, the 
statement of Theorem 3 is false. For example, the trivial solution y = to the differential 
equation y' = y is unstable on [xq, oo) for any xq > — oo. More generally, given the initial 
value problem 

y' = Mj , y{xo) = yo , 

with A a complex number, the solution y{x) = yoexp{X{x — xq)) is unstable for 3?A > 0; 
the solution is stable in the sense of Lyapunov for 3?A < and is asymptotically stable for 
3fiA < 0. 

In the next section we shall consider numerical methods for the approximate solution 
of the initial value problem (1-2). Since everything we shall say has a straightforward 
extension to the case of the system (5), for the sake of notational simplicity we shall restrict 
ourselves to considering a single ordinary differential equation corresponding to m = 1. We 
shall suppose throughout that the function / satisfies the conditions of Picard's Theorem 
on the rectangle R and that the initial value problem has a unique solution defined on the 
interval [xq,Xm], — oo < xq < Xm < oo. We begin by discussing one-step methods; this 
will be followed in subsequent sections by the study of linear multi-step methods. 

2 One-step methods 

2.1 Euler's method and its relatives: the ^-method 

The simplest example of a one-step method for the numerical solution of the initial value 
problem (1-2) is Euler's method^. 

Euler's method. Suppose that the initial value problem (1-2) is to be solved on the 
interval [xq, Xm]- We divide this interval by the mesh-points Xn = xo+nh, n = 0, . . . ,N, 
where h = (Xm — xq)/N and N is a positive integer. The positive real number h is called 
the step size. Now let us suppose that, for each n, we seek a numerical approximation y„ 
to y{xn), the value of the analytical solution at the mesh point Xn- Given that y{xo) = yo 

"Leonard Euler (1707-1783) 
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is known, let us suppose that we have aheady calculated y^, up to some n, < n < N — 1; 
we define 

Vn+l =yn + hfiXn, Vn) , n = 0, . . . , iV - 1 . 

Thus, taking in succession ?t- = 0, 1,...,A^— 1, one step at a time, the approximate values 
Un at the mesh points be easily obtained. This numerical method is known as 

Euler's method. 

A simple derivation of Euler's method proceeds by first integrating the differential 
equation (1) between two consecutive mesh points Xn and Xn+i to deduce that 

y{xn+i) = y{xn) + j f{x,y{x))dx, n = 0,...,N-l, (14) 

and then applying the numerical integration rule 

rXn+l 

/ g{x) dx hg{xn) , 
called the rectangle rule, with g{x) = f{x,y{x)), to get 

y{Xn+l) ^ y{x„) + hf{Xn,y{Xn)) , n = 0, ...iV-1, y{xo)=yo. 

This then motivates the definition of Euler's method. The idea can be generalised by 
replacing the rectangle rule in the derivation of Euler's method with a one-parameter 
family of integration rules of the form 

/■a; ,1+1 

/ g{x)dx^h[{l-e)g{xn) + eg{xn+i)] , (15) 

JXn 

with 6 G [0, 1] a parameter. On applying this in (14) with g{x) = f{x,y{x)) we find that 
y{Xn+l) ^ y{Xn)+h[{l-e)f{Xn,y{Xn))+Of{Xn+l,y{Xn+l))] , n = 0, . . . , AT - 1 , 

y{xo) = yo ■ 

This then motivates the introduction of the following one-parameter family of methods: 
given that yo is supplied by (2), define 

yn+i = yn + h[{l- 6)f{x„, yn) + Of{xn+\,yn+\)] , n = G,...,N -I . (16) 

parametrised by 6* G [0, 1]; (16) is called the ^-method. Now, for ^ = we recover Euler's 
method. For ^ = 1, and yo specified by (2), we get 

yn+l=yn + hf{Xn+uyn+l) , n = 0, . . . , iV - 1 , (17) 

referred to as the Implicit Euler Method since, unlike Euler's method considered above, 
(17) requires the solution of an implicit equation in order to determine yn+i, given 
In order to emphasize this difference, Euler's method is sometimes termed the Explicit 
Euler Method. The scheme which results for the value of ^ = 1/2 is also of interest: yo 
is supplied by (2) and subsequent values yn+i are computed from 

yn+l =yn+\h [/(Xn, yn) + /(Xn+1, yn+l)] , 71 = 0, ...,iV-l; 
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Vk iox = {) 


Vk for e = 1/2 


yk for 6* = 1 

















1 


0.1 





0.00500 


0.00999 


2 


0.2 


0.01000 


0.01998 


0.02990 


3 


0.3 


0.02999 


0.04486 


0.05955 


4 


0.4 


0.05990 


0.07944 


0.09857 



Table 1: The values of the numerical solution at the mesh points 



this is called the Trapezium Rule Method. 

The ^-method is an explicit method for ^ = and is an implicit method for < < 1, 
given that in the latter case (16) requires the solution of an implicit equation for yn+i- 
Further, for each value of the parameter 9 G [0, 1], (16) is a one-step method in the sense 
that to compute yn+i we only use one previous value Methods which require more 
than one previously computed value are referred to as multi-step methods, and will be 
discussed later on in the notes. 

In order to assess the accuracy of the ^-method for various values of the parameter 9 
in [0, 1], we perform a numerical experiment on a simple model problem. 

Example 1 Given the initial value problem y' = x — y"^ , y{0) = 0, on the interval of 
X G [0,0.4], we compute an approximate solution using the 9-method, for 9 = 0,9 = 1/2 
and 9 = 1, using the step size h = 0.1. The results are shown in Table 1. In the case 
of the two implicit methods, corresponding to 9 = 1/2 and 9 = 1, the nonlinear equations 
have been solved by a fixed-point iteration. 

For comparison, we also compute the value of the analytical solution y{x) at the mesh 
points Xn = 0.1 * n, n = 0, . . . ,4. Since the solution is not available in closed form,^ we 
use a Picard iteration to calculate an accurate approximation to the analytical solution on 
the interval [0, 0.4] and call this the "exact solution". Thus, we consider 



yo{x) = 



yk{x) 



yl- 



i(e)w^, fc = i,2,, 



Hence, 



yoix) = , 



yi{x) 
y2{x) 



-X 



1 



7:X^ — —X" 
2 20 



■'Using MAPLE V, wc obtain the solution in terms of Bessel functions: 
> dsolve({diff (y(x) ,x) + y(x)*y(x) = x, y(0)=0}, y(x)); 



■s/x 



'^V3BesselK(^, ^x3/2) _^ „ 
^— ^ Bessell( — , - x^/"^) 

TT 3 3 



y{x) = -- 



\ 



V3BesselK(J, ^x^/^) . o 
2-^ + Bessell(-, - xV^) 

TT 3 3 
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k 


Xk 


y{xk) 











1 


0.1 


0.00500 


2 


0.2 


0.01998 


3 


0.3 


0.04488 


4 


0.4 


0.07949 



Table 2: Values of the "exact solution" at the mesh points 



y^^^^ = 2^' - 20^' + 160^' - 4^^" • 
It is easy to prove by induction that 

y(^) = l^' - + - 7^^'' + o (x'^) , 

' 1 20 160 4400 V ; ' 

Tabulating ysix) on the interval [0,0.4] with step size h = 0.1, we get the "exact solution" 
at the mesh points shown in Table 2. 

The "exact solution" is in good agreement with the results obtained with 9 = 1/2: the 
error is < 5* 10~^. For 9 = and 9 = 1 the discrepancy between y^ and y{xk) is larger: it 
is < 3 * 10~^. We note in conclusion that a plot of the analytical solution can be obtained, 
for example, by using the MAPLE V package by typing in the following at the command 
line: 



> with(DEtools) : 

> DEplotCdiff (y(x) ,x)+y(x)*y(x)=x, y(x) , x=0..0.4, [[y(0)=0]], 
y=-0.1..0.1, stepsize=0.05) ; 

So, why is the gap between the analytical solution and its numerical approximation in 
this example so much larger for 9^1/2 than for 9 = 1/2? The answer to this question is 
the subject of the next section. 



2.2 Error analysis of the ^-method 

First we have to explain what we mean by error. The exact solution of the initial value 
problem (1 2) is a function of a continuously varying argument x G [xqi^m]) while the 
numerical solution y„ is only defined at the mesh points a:„, n = 0, . . . , A^, so it is a function 
of a "discrete" argument. We can compare these two functions either by extending in some 
fashion the approximate solution from the mesh points to the whole of the interval [xq, Xm\ 
(say by interpolating between the values or by restricting the function y to the mesh 
points and comparing y{xn) with y„ for n = 0, . . . , N . Since the first of these approaches 
is somewhat arbitrary because it does not correspond to any procedure performed in a 
practical computation, we adopt the second approach, and we define the global error e 

by 

en = y{xn) -yn , n = 0,...,N . 

We wish to investigate the decay of the global error for the ^-method with respect to the 
reduction of the mesh size h. We shall show in detail how this is done in the case of Euler's 
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method {9 = 0) and then quote the corresponding result in the general case (0 < < 1) 
leaving it to the reader to fill the gap. 

So let us consider Euler's explicit method: 

Vn+i = yn + hf{xn, Vn) , n = 0, . . . , iV , yo = given . 

The quantity 

Tn = ^(^-+^)^- ^(^--) _ y(^„)) , (18) 

obtained by inserting the analytical solution y{x) into the numerical method and dividing 
by the mesh size is referred to as the truncation error of Euler's explicit method and 
will play a key role in the analysis. Indeed, it measures the extent to which the analytical 
solution fails to satisfy the difference equation for Euler's method. 

By noting that f{xn-,y{xn)) = y'{xn) and applying Taylor's Theorem, it follows from 
(18) that there exists ^„ G {xn,Xn+i) such that 

Tn = 2 VCCn) , (19) 

where we have assumed that that / is a sufficiently smooth function of two variables so as 
to ensure that y" exists and is bounded on the interval [xq, Xm]- Since from the definition 
of Euler's method 

^ — 1^ J \Xn, yn) ) 

on subtracting this from (18), we deduce that 

+ h[f{Xn, y{Xn)) - f{Xn, Vn)] + hTn . 

Thus, assuming that \yn — yo\ ^ from the Lipschitz condition (3) we get 

|e„+i| < {l + hL)\en\+h\Tn\ , n = 0, . . . , iV - 1 . 
Now, let T = maxo<n<Ar-i \Tn\ ', then, 

|en+i| < {I + hL)\en\ + hT , n = 0,...,N -I . 
By induction, and noting that 1 + hL < e^^ , 

\en\ < ^[(l + /iL)"-l] + (l + /iL)»|eo| 

- I (e^^''""''°^ - l) + e-^(^"-^°) |eo| , n = 1, . . . , iV . 
This estimate, together with the bound 

\T\<hM2, M2= max ly"(x)| , 



which follows from (19), yields 



e„| < e-^(^"-^o)|eo| + ^ (e-^^^""^"^ - l) , n = 0,...,N. (20) 
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To conclude, we note that pursuing an analogous argument it is possible to prove that, 
in the general case of the ^-method, 



^11 It •^IT' \ 



h 



M2 + \hM3 



(21) 



for n = 0, . . . ,N, where now M3 = majz^^^^^xM] Iv (^) I- absence of rounding errors 

in the imposition of the initial condition (2) we can suppose that eo = yixo) — yo = 0. 
Assuming that this is the case, we see from (21) that |e„| = 0{h'^) for 9 = 1/2, while for 
6 = and 9 = 1, and indeed for any 9 7^ 1/2, |e„| = 0(h) only. This explains why in 
Tables 1 and 2 the values y„ of the numerical solution computed with the trapezium-rule 
method {9 = 1/2) were considerably closer to the analytical solution y{xn) at the mesh 
points than those which were obtained with the explicit and the implicit Euler methods 
{9 = and 9 = 1, respectively). 

In particular, we see from this analysis, that each time the mesh size h is halved, the 
truncation error and the global error are reduced by a factor of 2 when 9 7^ 1/2, and by a 
factor of 4 when 9 = 1/2. 

While the trapezium rule method leads to more accurate approximations then the 
forward Euler method, it is less convenient from the computational point of view given 
that it requires the solution of implicit equations at each mesh point Xn+i to compute 
yn+i- An attractive compromise is to use the forward Euler method to compute an initial 
crude approximation to y{xn+i) and then use this value within the trapezium rule to 
obtain a more accurate approximation for y{xn+i)- the resulting numerical method is 

yn+l=yn + ^h[f{Xn,yn)+f{Xn+l,yn + hf{Xn,yn))], U = 0, . . . , N , yo = given , 

and is frequently referred to as the improved Euler method. Clearly, it is an explicit 
one-step scheme, albeit of a more complicated form than the explicit Euler method. In the 
next section, we shall take this idea further and consider a very general class of explicit 
one-step methods. 

2.3 General explicit one-step method 

A general explicit one-step method may be written in the form: 

yn+i = yn + h^Xn, yn,h) , n = 0,...,N -1 , yo = y{xo) [= specified by (2)] , (22) 

where <!>(•, •; •) is a continuous function of its variables. For example, in the case of Euler's 
method, $(a;„, y„; h) = f{xn, yn), while for the improved Euler method 

$(Xn, yn, h) = ^ [fiXn, yn) + /(a^n + ^, 2/n + hf{Xn, yn))] ■ 

In order to assess the accuracy of the numerical method (22), we define the global 
error, e^, by 

~ y{,Xn) yn • 
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We define the truncation error, r„, of the method by 

Tn = ^(^"+^)^- y^'^-) _ . (23) 

The next theorem provides a bound on the global error in terms of the truncation 
error. 

Theorem 4 Consider the general one-step method (22) where, in addition to being a 
continuous function of its arguments, $ is assumed to satisfy a Lipschitz condition with 

respect to its second argument; namely, there exists a positive constant L$ such that, for 
< h < ho and for the same region R as in Picard 's theorem, 

\^{x,y;h)-^{x,z;h)\<L^\y-z\, for {x,y), {x, z) in R . (24) 

Then, assuming that \yn — yo\ < Ym, it follows that 



\e„\ < e^*(^"-^o)|eo| + 



n = 0,...,iV, (25) 



where T = maxo<n<Ar-i \Tn\ ■ 

Proof: Subtracting (22) from (23) we obtain: 

+ h[^{xn, y{xn); h) - $(a;„, h)] + hTn ■ 
Then, since {xn,y{xn)) and {xn,yn) belong to R, the Lipschitz condition (24) implies that 

|en+i| < |e„| + /iL$|e„| + h\Tn\ , n = 0, . . . ,N - 1 . 



That is. 
Hence 



|en+i| < (1 + hL^)\en\ + h\Tn\ , n = 0, . . . ,N - I . 

|ei| < {1 + hL^)\eo\ + hT , 

|e2| < {1 + hL^f\eo\ + h[l + {1 + hL^)]T , 

leal < {1 + hL^f\eo\ + h[l + {1 + hL^) + {1 + hL^f]T , 
etc. 

|e„| < (l + /iL$)"|eo| + [(l + /iL$)"-l]r/L$ . 

Observing that 1 + HL^ < exp(/iL$), wc obtain (25). o 

Let us note that the error bound (20) for Euler's explicit method is a special case 
of (25). We highlight the practical relevance of the error bound (25) by focusing on a 
particular example. 

Example 2 Consider the initial value problem y' = tan^-*^ y, y(0) = yo, and suppose that 
this is solved by the explicit Euler method. The aim of the exercise is to apply (25) to 
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quantify the size of the associated global error; thus, we need to find L and M2. Here 
f{x,y) = i&n~^y, so by the Mean Value Theorem 



\f{x,y) - f{x,z)\ 



where 77 lies between y and z. In our case 

df 



dy 



and therefore L = 1. To find M2 we need to obtain a bound on \y"\ (without actually 
solving the initial value problem!). This is easily achieved by differentiating both sides of 
the differential equation with respect to x: 

y" = ^(tan-i y) = (1 + y^'^ = (1 + vY' tan-^ y . 

Therefore \y"{x)\ < M2 = Inserting the values of L and M2 into (20), 

|en| < e^"|eo| + ^7r(e''" - 1) /i, n = 0,...,N . 

In particular if we assume that no error has been committed initially (i.e. eo = 0), we 
have that 

|e„| < -7r(e^" - n = (},...,N. 

Thus, given a tolerance TOL specified beforehand, we can ensure that the error between 
the (unknown) analytical solution and its numerical approximation does not exceed this 
tolerance by choosing a positive step size h such that 

h < -(e^^ - ly^TOL . 

TT 

For such h we shall have \y{xn) — yn\ = \&n\ ^ TOL for each n = 0, . . . , N , as required. 
Thus, at least in principle, we can calculate the numerical solution to arbitrarily high ac- 
curacy by choosing a sufficiently small step size. In practice, because digital computers use 
finite-precision arithmetic, there will always be small (but not infinitely sm,all) pollution ef- 
fects due to rounding errors; however, these can also be bounded by performing an analysis 
similar to the one above where /(x^, yn) is replaced by its finite-precision representation. 

Returning to the general one-step method (22), we consider the choiee of the function 
$. Theorem 4 suggests that if the truncation error 'approaches zero' as h ^ then the 
global error 'converges to zero' also (as long as |eo| — > when h 0). This observation 
motivates the following definition. 



Definition 2 The numerical method (22) is consistent with the differential equation 
(1) if the truncation error defined by (23) is such that for any e > there exists a positive 
h{e) for which |T„| < e for < h < h{e) and any pair of points {xn, y{xn)), {xn+i, y{xn+i)) 
on any solution curve in R. 
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For the general one-step method (22) we have assumed that the function $(•,•;•) is 
continuous; also y' is a continuous function on [xc-'^m]- Therefore, from (23), 

lim Tn = y'{Xn) - ^{Xn, y{Xn)] 0) . 

h—*0 

This implies that the one-step method (22) is consistent if and only if 

^x,y;0) = f{x,y). (26) 

Now we are ready to state a convergence theorem for the general one-step method (22). 

Theorem 5 Suppose that the solution of the initial value problem (1-2) lies in R as does 
its approximation generated from (22) when h < ho. Suppose also that the function $(•, •; •) 
is uniformly continuous on R x [0,ho] and satisfies the consistency condition (26) and the 
Lipschitz condition 



\^{x, y; h) — $(x, z\h)\ < L^\y — z\ on R x [0, ho] ■ 



(27) 



Then, if successive approximation sequences (y„), generated for Xn = xq + nh, n = 
1,2, . . . , N , are obtained from (22) with successively smaller values of h, each less than Hq, 
we have convergence of the numerical solution to the solution of the initial value problem 
in the sense that 







as h ^ 0, Xn X & [xq, Xm] ■ 



Proof: Suppose that = {Xm — xo)/N where A'^ is a positive integer. We shall assume 
that N is sufficiently large so that h < hg. Since y{xQ) = yo and therefore cq = 0, Theorem 
4 implies that 



\y{xn) -yn\< 



^L^{Xm—xi.>) _ 1 
La, 



max 

0<m<n-l 



n = l,...,7V. 



(28) 



Prom the consistency condition (26) we have 

'y{xn+i) - y{xn) 



h 



f{Xn,y{Xn)) 



+ [^{Xn, y{Xn);0) - ^{Xn, y{Xn)\ h)] . 



According to the Mean Value Theorem the expression in the first bracket is equal to 

y'iO - y'ixn), where ^ G Since y'(-) = f{-,y{-)) = $(-,y(-);0) and $(•,•:•) is 

uniformly continuous on R x [0, h^], it follows that y' is uniformly continuous on [xq, Xm]- 
Thus, for each e > there exists /ii(e) such that 



\y'iO - y'ixn)\ < for h < hlie), n = 0, 1, . 



Also, by the uniform continuity of $ with respect to its third argument, there exists /12(e) 
such that 



MXn, y{Xny,0) - HXn, y{Xn); h)\ < -6 for k < /12(e), n = 0, 1, . . . , N - 1 . 
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Thus, defining h{e) = min(/ti(e), /12(e)), we have 



Tn\<€ 



ioT h < h{e), n = 0,l,...,N -1 . 



Inserting this into (28) we deduce that \y{xn) — yn\ as h ^ 

\y{x) - yn\ < \y{x) - y{xn)\ + \yixn) - y-, 



In 



0. Since 



and the first term on the right also converges to zero as /i — >■ by the uniform continuity 
of y on the interval [xc^m], the proof is complete, o 

We saw earlier that for Euler's method the absolute value of the truncation error T„ 
is bounded above by a constant multiple of the step size h, that is 



where is a positive constant, independent of h. However there are other one-step 
methods (a class of which, called Runge-Kutta methods, will be considered below) for 
which we can do better. More generally, in order to quantify the asymptotic rate of decay 
of the truncation error as the step size h converges to zero, we introduce the following 
definition. 

Definition 3 The numerical method (22) is said to have order of accuracy p, if p is 
the largest positive integer such that, for any sufficiently smooth solution curve {x,y{x)) 
in R of the initial value problem (1-2), there exist constants K and ho such that 

\Tn\<KhP forO<h<ho 

for any pair of points {xn,y{xn)), (.x„_|_i, y(x„+i)) on the solution curve. 

Having introduced the general class of explicit one-step methods and the associated 
concepts of consistency and order of accuracy, we now focus on a specific family: explicit 
Runge-Kutta methods. 

2.4 Runge-Kutta methods 

In the sense of Definition 3 Euler's method is only first-order accurate; nevertheless, it 
is simple and cheap to implement because to obtain yn+i from y„ we only require a 
single evaluation of the function f at (x^, ?/„). Runge-Kutta methods aim to achieve 
higher accuracy by sacrificing the efficiency of Euler's method through re-evaluating /(•, •) 
at points intermediate between (x„,y(x„)) and {xn+i,y{xn+i))- The general i?-stage 
Runge-Kutta family is defined by 



|T„| < Kh 



ior <h<ho 



Vn+l 



yn + h^{xn,yn;h) , 



Hx,y;h) 




f{x,y) , 



'r 




(29) 



r-l 



r = 2,.... 



R . 
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a = Be 


B 




T 




c 



where e = (1, . . . , l)"*^. 
Figure 1: Butcher table of a Runge-Kutta method 



In compressed form, this information is usually displayed in the so-called Butcher table 

displayed in Figure 1. 

One-stage Runge— Kutta methods. Suppose that R = 1. Then, the resulting one- 
stage Runge-Kutta method is simply Euler's explicit method: 

Vn+l =yn + hf{Xn, Vn) ■ (30) 

Two-stage Runge Kutta methods. Next, consider the case of i? = 2, corresponding 
to the following family of methods: 

Vn+l =yn + h{ciki + C2k2) , (31) 

where 

^1 = f{Xn,yn) , (32) 

k2 = fixn + a2h, Un + b2ihki) , (33) 

and where the parameters ci, C2, 02 and 621 are to be determined.^ Clearly (31-33) can 
be rewritten in the form (22) and therefore it is a family of one step methods. By the 
condition (26) , a method from this family will be consistent if and only if 

Ci + C2 = 1 . 

Further conditions on the parameters are obtained by attempting to maximise the order 
of accuracy of the method. Indeed, expanding the truncation error of (31-33) in powers 
of h, after some algebra we obtain 



r„ = lhy"iXn) + lhY'iXn) 

-C2h[a2fx + 621/y/] - C2/t^ 



^alfxx + a2b2ifxyf + 2^21/2/2//^ 



Here we have used the abbreviations / = /(a;„, y{xn)), fx = '^{xn, y{xn)), etc. On noting 
that y" = fx + fyf, it follows that T„ = Oiyh?) for any / provided that 

a2C2 = b2lC2 = ^ , 

which implies that if 621 = 02, C2 = l/(2a2) and ci = 1 — l/(2a2) then the method is 
second-order accurate; while this still leaves one free parameter, 02, it is easy to see that 
no choice of the parameters will make the method generally third-order accurate. There 
are two well-known examples of second-order Runge-Kutta methods of the form (31-33): 



^We note in passing that Euler's method is a member of this family of methods, corresponding to ci = 1 
and C2 = 0. However we are now seeking methods that are at least second-order accurate. 
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a) The modified Euler method: In this case we take 02 = ^ to obtain 

Vn+l =yn + hf (xn + ^h, Vn + ^hf{Xn, Vnij ; 

b) The improved Euler method: This is arrived at by choosing 02 = 1 which gives 

+ 2^ [fi^n, Vn) + f{Xn + h,yn + hf{Xn, Vn))] ■ 

For these two methods it is easily verified by Taylor series expansion that the truncation 
error is of the form, respectively, 



Tn = ^h' 
6 



fyFi + 
fyFi - \f2 



where 



Fl= ffy and F2 = Ux + 2fUy + f fyy 



The family (31-33) is referred to as the class of explicit two-stage Runge-Kutta methods. 

Exercise 1 Let a he a non-zero real number and let x„ = a + nh, n = 0, . . . , N , be a 
uniform mesh on the interval [a, b] of step size h = {b — a)/N. Consider the explicit one- 
step method for the numerical solution of the initial value problem y' = f{x, y), y(a) = yo, 
which determines approximations yn to the values y{xn) from the recurrence relation 

Vn+l =yn + h{l - a)f{Xn,yn) + haf (xn + -^,yn + ^f{Xn,yn, 



2a 



2a' 



Show that this method is consistent and that its truncation error, T„{h,a), can be 
expressed as 



Tnih,a) 



80; 



a-l] y"'{Xn) + y"iXn)-^iXn,y{Xn)) 



+ 0(/i^) 



This numerical method is applied to the initial value problem y' = —y^, y(0) = 1, where 

p is a positive integer. Show that if p = 1 then Tn{h,a) = Oiji^) for every non-zero real 
number a. Show also that if p > 2 then there exists a non-zero real number ao such that 
Tn{h,ao) = 0{h^). 



Solution: Let us define 



Hx,y;h) = {l-a)f{x,y) + af [ x + ^. jj + ^f{x,y) ] . 



2a 

Tlien the numerical method can be rewritten as 

Vn+l =yn + h^{Xn,yn\h) 

Since 

$(a;,j/;0) = f{x,y) , 



2a 



15 



the method is consistent. By definition, the truncation error is 

Tn[h, a) = ^{Xn, y{Xn); h) . 

n 

Wc! shall perform a Taylor expansion of T„(/i, a) to show that it can be expressed in the desired 
form. Indeed, 

h 

Tn{h,a) = y'{xn) + -y"ixn) + —y"'{xn) 

-(1 - a)y'{xn) - af{xn + ^^yi^n) + ^y'M) + 0{h^) 

h h? 
= y'{xn) + ^y"{xn) + yy"'(a;„) - (1 - a)y'{xn) 



-a 

a 
~2 



f{Xn, y{Xn)) + ^f^i^n, y{Xn)) + ^fvi^n, y{Xn))y' {Xn) 



h_ 
2a 



2 / , \ 2 



fxx{Xn,y{Xn)) + 2 7— fxy{Xn,y{Xn))y' {Xn) 



2a 



+ 



_h_ 
2a 



fyy{Xn,y{Xn))[y'{Xn)]'^ 



+ 0{h^) 



y'{xn) - (1 - a)y'{xn) - ay'{xn) 

+ '^y"iXn) - ^ [fx{Xn,y{Xn)) + fy{Xn, y{Xn))y' (Xn)] 

+ -^y"'{Xn) - — [fxx{Xn,y{Xn)) + 2fxv{Xn,y{Xn))y'{Xn) 

+ fyy{Xn,y{Xn))[y'iXn)?] + 0{h^) 
'^y"'{Xn) - ^ [y"'{Xn) - y"{Xr,)fy{Xr„y{Xn))]+0{h^) 



8a 



-a - 1 ) 2/"'(a;„) + y"(a;„) — (a;„,y(x„)) 



as required. 

Now let us apply the method to y' = —y^, with p>l. If p = 1, then y'" = —y" = y' = —y, so 
that 

/,2 

T^{h,a) = -—y{xn)+0{h'') . 



As y{xn) = e ^ 0, it follows that 



for all (non-zero) a. 

Finally, suppose that p>2. Then 



Tn{h,a) = 0{h'') 



and 

and therefore 



Tn{h, a) = - 



y'" = p{2p - l)y^P-^y' = -p{2p - , 

y^P-\xn) + 0{h^). 



8a 



a-l]p{2p-l)+p^ 
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Choosing a such that 
namely 



\q.-\\ p(2p- 1) = 



3p- 3 

OL = aQ = , 

8p — 4 

gives 

T„(/i,ao)=0(/i3). 

We note in passing that for p > 1 the exact solution of the initial value problem y' = — y^, 
2/(0) = 1, is = \{v - \)x + l^l^^-^y o 
Three-stage Runge-Kutta methods. Let us now suppose that = 3 to illustrate 
the general idea. Thus, we consider the family of methods: 

=yn + h [cifci + 02^2 + c^k^] , 

where 

^1 = f{x,y), 

k2 = fix + ha2,y + hb2iki), 

h = fix + ha3, y + ^631^1 + ^632^2) , 

02 = 621 , 03 = ^31 + &32 • 

Writing 621 = 02 and 631 = 03 — 632 in the definitions of k2 and respectively and 
expanding k2 and ^3 into Taylor series about the point (x, y) yields: 

k2 = f + ha2ifx + kify) + ^h^al{f^^ + 2kiUy + kjfyy) + 0{h^) 

= f + ha^iU + ffy) + \h^al{U, + 2fUy + f fyy) + 0(/i3) 
= / + ha2Fi + ^h?alF2 + 0{h^) , 



where 
and 



Fl = fx + ffy and F2 = fxx + 2ffxy + f^fyy , 



ks = f + h {as fx + [(as - 532)/ci + 632A;2] fy} 

+ ^^^ {(4fxx + 2a3 [(03 - 632)A;i + 632^2] fxy 

+ [{as - b32)kl + 632^:2]' fyy} + 0{h^) 

= / + /laaFi + (a2h2Fify + \alF^ + 0(/i') . 

Substituting these expressions for k2 and A;3 into (29) with i? = 3 we find that 

^{x,y,h) = (ci + C2 + C3)/ + /i(c2a2 + 0303)^1 

+^h^ [20302632^1/2, + {c2al + C3ai) F2] + 0{h^) . (34) 
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We match this with the Taylor series expansion: 

y^'^ + ^l-y^'^^ = y'{x) + lhy"ix) + lhV'{x) + 0{h') 
= f+^hF, + ^h''{F^fy + F2) + 0{h^). 

This yields: 

C1 + C2 + C3 = 1 , 
1 

C2a2 + 0303 = - , 
0202 + C30I = ^ , 

0302^32 = ^ • 
6 

Solving this system of four equations for the six unknowns: ci, C2, C3, 02, 03, 632, we obtain 
a two-parameter family of 3-stage Runge-Kutta methods. We shall only highlight two 
notable examples from this family: 

(i) Heun's method corresponds to 

1 3 12 2 

ci = ^ , C2 = , C3 = - , «2 = 2 ' "3 = 2, ^'32 = 2 ' 

yielding 

Vn+i = yn + ^h{ki + 3k3) , 

= f{Xni Un) 5 

h = f (xn + ^h, yn + ^hki 

/ 2 2 

h = f {xn + -h, yn + -hk2 

(ii) Standard third-order Runge— Kutta method. This is arrived at by selecting 

12 11 

Cl = g , C2 = - , C3 = - , ^2 = 2 ' «3 = 1 , O32 = 2 , 



yielding 



Vn+l = Vn + ^h{ki+ 4k2 + ks) , 



kl = f{Xn,yn) 

k2 = f (xn + ^h, yn + ^hki ) , 
h = f {xn + h,yn- hki + 2hk2) 
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Four-stage Runge— Kutta methods. For i? = 4, an analogous argument leads to a 
two-parameter family of four-stage Runge-Kutta methods of order four. A particularly 
popular example from this family is: 



Vn+i = yn + ^h {ki + 2k2 + 2k3 + ki) , 

where 

^1 = fi^ni Un) ) 

k2 = f (^Xn + ^h, yn + ^hki 

kz = J {xn + ^h, yn + ^hk2 
ki = f{xn + h,yn + hk3) . 

Here k2 and ^3 represent approximations to the derivative y'{-) at points on the solution 
curve, intermediate between {xn, y{xn)) and (xn+i,y{xn+i)), and $(x„, yn', h) is a weighted 
average of the ki, i = 1, . . . , 4, the weights corresponding to those of Simpson's rule method 
(to which the fourth-order Runge-Kutta method reduces when ^ = 0). 

In this section, we have constructed i?-stage Runge-Kutta methods of order of accuracy 
0{h^), R = 1,2, 3, 4. Is is natural to ask whether there exists an R stage method of order 
R for R> 5. The answer to this question is negative: in a series of papers John Butcher 
showed that for R = 5,6, 7, 8, 9, the highest order that can be attained by an i?-stage 
Runge-Kutta method is, respectively, 4, 5, 6, 6, 7, and that for i? > 10 the highest order is 
< R-2. 



2.5 Absolute stability of Runge-Kutta methods 

It is instructive to consider the model problem 

y' = Xy, y(0) = yo (7^ 0), (35) 

with A real and negative. Trivially, the analytical solution to this initial value problem, 
y{x) = ?/oexp(Ax), converges to at an exponential rate as x — > -|-oo. The question that 
we wish to investigate here is under what conditions on the step size h does a Runge-Kutta 
method reproduce this behaviour. The understanding of this matter will provide useful 
information about the adequate selection of h in the numerical approximation of an initial 
value problem by a Runge-Kutta method over an interval [xq, Xm] with Xm » xq. For 
the sake of simplicity, we shall restrict our attention to the case of i?-stage methods of 
order of accuracy R, with 1 < < 4. 

Let us begin with R = 1. The only explicit one-stage first-order accurate Runge-Kutta 
method is Euler's explicit method. Applying (30-35) yields: 

yn+i = (1 + h)yn , n>o , 

where h = Xh. Thus, 

yn = (1 + hTyo . 
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Consequently, the sequence {?/n}^o ^^^^ converge to if and only if |1 + ^| < 1, yielding 
h G (—2,0); for such h the Eulcr's explicit method is said to be absolutely stable and 
the interval (—2, 0) is referred to as the interval of absolute stability of the method. 
Now consider R = 2 corresponding to two-stage second-order Runge-Kutta methods: 

Vn+l =yn + h{ciki + C2k2) , 

where 

ki = f{xn,yn), k2 = f{xn + a2h,yn + b2ihki) 
with ^ 

Ci+C2 = l , a2C2 = 621C2 = 2 • 

Applying this to (35) yields, 

Vn+i = (l + h+ ^hA yn , n > , 



and therefore 



\jn= (l + h+ ^hA yo 



Hence the method is absolutely stable if and only if 

<1, 



1 + h+lh'' 



namely when h G (—2, 0). 

In the case of i? = 3 an analogous argument shows that 

yn+i= (l + h+^h^ + ^hAyn . 



Demanding that 



< 1 



then yields the interval of absolute stability: h G (—2.51,0). 
When R = 4, we have that 

Vn+l =(^l + h+^P + ^h^ + ^h^^ yn , 

and the associated interval of absolute stability is ^ G (—2.78,0). 

For i? > 5 on applying the Runge-Kutta method to the model problem (35) still results 
in a recursion of the form 

Vn+l = AR{h)yn , n > , 

however, unlike the case when R = 1,2,3,4, in addition to h the expression Aji(h) also 
depends on the coefficients of the Runge-Kutta method; by a convenient choice of the free 
parameters the associated interval of absolute stability may be maximised. For further 
results in this direction, the reader is referred to the book of J.D. Lambert. 
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3 Linear multi-step methods 



While Runge-Kutta methods present an improvement over Euler's metliod in terms of 
accuracy, this is achieved by investing additional computational effort; in fact, Runge- 
Kutta methods require more evaluations of /(•, •) than would seem necessary. For example, 
the fourth-order method involves four function evaluations per step. For comparison, by 
considering three consecutive points Xn = Xn-i+h, Xn+i = Xn-i+2h, integrating the 
differential equation between Xn-i and Xn+i, and applying Simpson's rule to approximate 
the resulting integral yields 

y{xn+i) = y{xn-i) + / f{x,y{x))dx 

~ yiXn-l) + [f{Xn-l, y{Xn-l)) + 4/(Xn, y{Xn)) + 

which leads to the method 

2/n+l = yn-l + [f{Xn-l, J/n-l) + 4/(x„, ?/n) + /(^n+l, Z/n+l)] ■ (36) 

In contrast with the one-step methods considered in the previous section where only a 
single value was required to compute the next approximation ?/„+i, here we need two 
preceding values, y„ and yn-i to be able to calculate yn+i, and therefore (36) is not a 
one-step method. 

In this section we consider a class of methods of the type (36) for the numerical solution 
of the initial value problem (1-2), called linear multi-step methods. 

Given a sequence of equally spaced mesh points (xn) with step size h, we consider the 
general linear A;-step method 

k k 

Yl "iV^+j = ^ H Pjfi^n+j, yn+j) , (37) 
j=0 j=0 

where the coefficients aQ,...,ak and /3o,---,/?A; are real constants. In order to avoid 
degenerate cases, we shall assume that ak ^ and that ao and /3o are not both equal to 
zero. If /?fc = then yn+k is obtained explicitly from previous values of yj and f{xj,yj), 
and the A;-step method is then said to be explicit. On the other hand, if ^ then yn+k 
appears not only on the left-hand side but also on the right, within f{xn+k,yn+k)) due 
to this implicit dependence on yn+k the method is then called implicit. The numerical 
method (37) is called linear because it involves only linear combinations of the {y„} and 
the {f{xn,yn)}', for the sake of notational simplicity, henceforth we shall write /„ instead 

of f{Xn,yn)- 

Example 3 We have already seen an example of a linear 2-step method in (36); here we 
present further examples of linear multi-step methods. 

a) Euler's method is a trivial case: it is an explicit linear one-step method. The im- 
plicit Euler method 

yn+l =yn + hf{Xn+l,yn+l) 

is an implicit linear one-step method. 
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b) The trapezium method, given by 

Vn+l =yn + ^h[fn+l + fn] 

is also an implicit linear one-step method. 

c) The four-step Adams^- Bashforth method 

yn+4 = yn+3 + ^/l[55/„+3 - 59/„+2 + 37/„+i - 9/„] 

is an example of an explicit linear four-step method; the four-step Adams - Moul- 
ton method 

yn+4: = yn+3 + ^^[9/n+4 + 19/n+3 - 5/„+2 - 9/n+l] 

is an implicit linear four-step method. 

The construction of general classes of linear multi-step methods, such as the (implicit) 
Adams- Bashforth family and the (explicit) Adams-Moulton family will be discussed in the 
next section. 

3.1 Construction of linear multi-step methods 

Let us suppose that tt„, n = 0, 1, . . ., is a sequence of real numbers. We introduce the shift 
operator E, the forward difference operator A+ and the backward difference operator A_ 

by 

E : Un^ Un+l , A+ : Un 1-^ (Un+l - Uji) , A- : Un {Un - Un-l) ■ 

Further, we note that E~^ exists and is given by E"'^ : Un+i i— »■ «n- Since 

A+ = E-I = EA_, A^=I-E-^ and £; = (/-A_)"\ 
where I signifies the identity operator, it fohows that, for any positive integer k, 



and 



Alun = {E- Ifun = f ^ ) U^+k 

j=0 V / 

Atun ={I- E-^fun = (^.]un- 

.7=0 V / 



3 



Now suppose that n is a real-valued function defined on R whose derivative exists and 
is integrable on [xo,a;„] for each n > 0, and let Un denote u{xn) where Xn = xq + nh, 
n = 0,1,..., are equally spaced points on the real line. Letting D denote d/dx, by 
applying a Taylor series expansion we find that 

{E^u)n = u{Xn + sh) = Un + sh{Du)n + ^{skf {D'^u)n + ■ ■ ■ 



oo 1 

k=o 



■^J. C. Adams (1819-1892) 
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and hence 



Thus, formally, 

hD = \nE = -ln(/ - A_) 
and therefore, again by Taylor series expansion, 

hu'{xn) = f A_ + + + 



Now letting u[x) = y{x) where y is the solution of the initial-value problem (1—2) and 
noting that u'{x) = y'{x) = f{x,y{x)), we find that 

hf{xn,yixn)) = (^A_ + ^Al + ^Al + y{xn) . 
On successive truncation of the infinite series on the right, we find that 

y{Xn) - y{Xn-l) ~ hf{Xn,y{Xn)) , 

3 1 

^y{Xn) - 2y(x„_i) + -y{Xn-2) ~ hf(Xn,y{Xn)) , 

11 3 1 

-^yi^ri) - 3y(Xn_i) + -y{Xn-2) - -y(x„_3) hf{Xn, y{Xn)) , 

and so on. These approximate equalities give rise to a class of implicit linear multi-step 
methods called backward differentiation formulae, the simplest of which is Euler's 
implicit method. 
Similarly, 

E-\hD) = hDE-^ = {I- A_)(-ln(7 - A_)) = -(/ - A_)ln(/ - A_) , 
and therefore 

hu'ixn) = (^A_ - -Al - -At + . . Un+l . 

Letting, again, u{x) = y{x) where y is the solution of the initial-value problem (1-2) and 
noting that u'{x) = y'{x) = f{x,y(x)), on successive truncation of the infinite series on 
the right results in 

y{Xn+l) - y{Xn) ~ hf{Xn,yiXn)) , 
^yi^n+l) - -^y{Xn-l) ^ hf{Xn,y{Xn)), 

^y{Xn+l) + ^y{Xn) - y{Xn-l) + ^yiXn-2) ~ hf{Xn,y{Xn)) , 

and so on. The first of these yields Euler's explicit method, the second the so-called 
explicit mid-point rule, and so on. 

Next we derive additional identities which will allow us to construct further classes of 
linear multi-step methods. Let us define 



-1 f^" 
D u{xn) = u{xo) + / u{^) d^ , 

JXQ 
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and observe that 

{E - I)D-\{xn) = ^"^'n(Od^ 

Now, 

{E-I)D-^ = A+D-^ = EA-D-^ = hEA-{hD)-'^ 

= -hEA^[ln{I - A^)]-K (38) 

Furthermore, 

{E-I)D-^ = EA.D-^ = A_ED-^ = A_{DE-^)-^ = hA_{hDE-^)-^ 

= -/iA_ [(I- A_)ln(/- A_)]-i . (39) 

Letting, again, u{x) = y{x) where y is the solution of the initial-value problem (1-2), 
noting that u'{x) = y'{x) = /(x, y{x)) and using (38) and (39) we deduce that 

y{Xn+l) - V{Xn) = / y'{i)di={E-I)D-^y'{Xn) = {E-I)D-^f{Xn,y{Xn)) 

^ \ -hEA-[\n{I-A-)\-^f{XnM^n)) .4.. 

\ -hA_ [{I - A_)ln(/ - A_)]-i /(x„, y{xn)) • ^ ^ 

On expanding ln(7 — A_) into a Taylor series on the right-hand side of (40) we find that 



y{xn+i) - y{xn) ~ h 
and 

y(Xn+l) - y{Xn) ^ h 



/ A_ A'' A"* A^ 

2 12 - 24 - 720 " 



/+1a_ + Aa2 +-a3 + — a^ + 

2 12 " 8 " 720 " 



f{xn,y{xn)) (41) 



f{Xn,y{Xn)). (42) 



Successive truncation of (41) yields the family of Adams-Moulton methods, while similar 
successive truncation of (42) gives rise to the family of Adams-Bashforth methods. 

Next, we turn our attention to the analysis of linear multi-step methods and introduce 
the concepts of stability, consistency and convergence. 

3.2 Zero-stability 

As is clear from (37) we need k starting values, yo, . . . , yk-i, before we can apply a linear k- 
step method to the initial value problem (1-2): of these, yo is given by the initial condition 
(2), but the others, yi, ■ ■ ■ ,yk-i, have to be computed by other means: say, by using a 
suitable Runge-Kutta method. At any rate, the starting values will contain numerical 
errors and it is important to know how these will affect further approximations y^, n > k, 
which are calculated by means of (37). Thus, we wish to consider the 'stability' of the 
numerical method with respect to 'small perturbations' in the starting conditions. 

Definition 4 A linear k-step method for the ordinary differential equation y' = f{x,y) 
is said to be zero-stable if there exists a constant K such that, for any two sequences 
(yn) and {fjn) which have been generated by the same formulae hut different initial data 
yo, yi, ■ ■ ■ ,yk-i and yo, yi, . ■ ■ ,yk-i, respectively, we have 

\yn - yn\ < max{|yo - yo\, |yi - yi\, \yk-i - yk-i\} (43) 
for Xn < Xm, and as h tends to 0. 
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We shall prove later on that whether or not a method is zero-stable can be determined 
by merely considering its behaviour when applied to the trivial differential equation y' = 0, 
corresponding to (1) with f{x, y) = 0; it is for this reason that the kind of stability 
expressed in Definition 4 is called zero stability. While Definition 4 is expressive in the 
sense that it conforms with the intuitive notion of stability whereby "small perturbations 
at input give rise to small perturbations at output", it would be a very tedious exercise 
to verify the zero-stability of a linear multi-step method using Definition 4 only; thus we 
shall next formulate an algebraic equivalent of zero-stability, known as the root condition, 
which will simplify this task. Before doing so we introduce some notation. 

Given the linear fc-step method (37) we consider its first and second characteristic 
polynomial, respectively 

k 

j=0 

k 

a{z) = Y.(3,z\ 

j=0 

where, as before, we assume that 

Now we are ready to state the main result of this section. 

Theorem 6 A linear multi-step method is zero-stable for any ordinary differential equa- 
tion of the form (1) where f satisfies the Lipschitz condition (3), if and only if its first 
characteristic polynomial has zeros inside the closed unit disc, with any which lie on the 
unit circle being simple. 

The algebraic stability condition contained in this theorem, namely that the roots of 
the first characteristic polynomial lie in the closed unit disc and those on the unit circle 
are simple, is often called the root condition. 

Proof: Necessity. Consider the linear /c-step method, applied to y' = 0: 

c^kVn+k + afc-i2/n+fc-i + . • • + ai2/„+i + a^yn = . (44) 
The general solution of this fcth order linear difference equation has the form 

Vn = Y,Ps{n)z-^ , (45) 

s 

where Zg is a zero of the first characteristic polynomial p{z) and the polynomial Ps{ ) 
has degree one less than the multiplicity of the zero. Clearly, if |2;s| > 1 then there are 
starting values for which the corresponding solutions grow like l^sl" and if = 1 and 
its multiplicity is rus > 1 then there are solutions growing like n'^''~^. In either case there 
are solutions that grow unbounded as n — > oo, i.e. as h ^ with nh fixed. Considering 
starting data yo, yi, . . . ,yk-i which give rise to such an unbounded solution (y„), and 
starting data yo = yi = ■ ■ ■ = Vk-i = for which the corresponding solution of (44) is (y„) 
with i/n = for all n, we see that (43) cannot hold. To summarise, if the root condition 
is violated then the method is not zero-stable. 
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Sufficiency. The proof that the root condition is sufficient for zero-stabiHty is long and 

technical, and will be omitted here. For details, sec, for example, K.W. Morton, Numerical 
Solution of Ordinary Differential Equations, Oxford University Computing Laboratory, 
1987, or P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, 
New York, 1962. o 



Example 4 We shall consider the methods from Example 3. 

a) The explicit and implicit Euler methods have first characteristic polynomial p{z) = 
z — 1 with simple root z = 1, so both methods are zero-stable. The same is true of 
the trapezium method. 

b) The Adams-Bashforth and Adams-Moulton methods considered in Example 3 have 
the same first characteristic polynomial, p{z) = z^{z — l), and therefore both methods 
are zero-stable. 

c) The three-step (sixth order accurate) linear multi-step method 

ll2/n+3 + 27y„+2 - 27yn+i - Ilyn = 3/l[/n+3 + 9/n+2 + 9/„+i + /„] 

is not zero-stable. Indeed, the associated first characteristic polynomial p{z) = 112;^+ 
27z^ — 272; — 11 has roots at z\ = 1, 22 ~ —0.3189, 23 k, —3.1356, so l^al > 1. 



3.3 Consistency 

In this section we consider the accuracy of the linear fc-step method (37). For this pur- 
pose, as in the case of one-step methods, we introduce the notion of truncation error. 
Thus, suppose that y{x) is a solution of the ordinary differential equation (1). Then the 
truncation error of (37) is defined as follows: 

rj, ^ Et=o Ky(a^n+j) - h(3jy'{Xn+j)] 

Of course, the definition requires implicitly that <t(1) = Ej=o/?i 7^ 0. A gain, as in the 
case of one-step methods, the truncation error can be thought of as the residual that is 
obtained by inserting the solution of the differential equation into the formula (37) and 
scaling this residual appropriately (in this case dividing through by hJ2j=o Pj) so that T„ 
resembles y' — f{x, y{x)). 

Definition 5 The numerical scheme (37) is said to be consistent with the differential 
equation (1) if the truncation error defined by (46) is such that for any e > there exists 
h{e) for which 

\Tn\ <e forO<h< h{e) 

and any {k + I) points (x„, y{xn)), ■ ■ ■ , {xn+k, y{xn+k)) on any solution curve in R of the 
initial value problem (1-2). 
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Now let us suppose that the solution to the differential equation is sufficiently smooth, 

and let us expand y{xn+j) and y'{xn+j) into a Taylor series about the point Xn and 
substitute these expansions into the numerator in (46) to obtain 



r„. = - J^ [Coy{xn) + Cihy'ixn) + C2hy{xn) + • ■ • ] (47) 



where 



k 

" 3 ■> 



Co = E« 

3=0 

k k 

Ci = 

i=i i=o 

fc -2 k 

i=i i=i 

etc. 

For consistency we need that ^ {) as h ^ Q and this requires that Co = and Ci = 0; 
in terms of the characteristic polynomials this consistency requirement can be restated in 
compact form as 

p{l) = and p'{l) = a{l) / . 

Let us observe that, according to this condition, if a linear multi-step method is consistent 
then it has a simple root on the unit circle at z = 1; thus the root condition is not violated 
by this zero. 

Definition 6 The numerical method (37) is said to have order of accuracy p if p is 
the largest positive integer such that, for any sufficiently smooth solution curve in R of the 
initial value problem (1-2), there exist constants K and Hq such that 

\Tn\ < KhP forO<h<ho 

for any {k + 1) points {xn, y{xn)), ■ ■ ■ , {xn+k: y{xn+k)) on the solution curve. 

Thus we deduce from (47) that the method is of order of accuracy p if and only if 

Co = Ci = . . . = Cp = and Cp+i ^ . 

In this case, 



Cp+i 



hPy^P+^\xn) + 0{hP+^) ; 



" a(l) 

the number Cp+i (7^ 0) is called the error constant of the method. 

Exercise 2 Construct an implicit linear two-step method of maximum order, containing 
one free parameter. Determine the order and the error constant of the method. 
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Solution: Taking ao = a as parameter, the method has the form 

yn+2 + OtlVn+l + ayn = h{P2fn+2 + /?l/n+l + /^o/n) , 

with a2 = 1, ao = a, P2 7^ 0. We have to determine four unknowns: ai, (32, Po, so we require 
four equations; these will be arrived at by demanding that the constants 

Co = ao + ai + a2 , 

Ci = ai+2-(/3o + /3i+/32) , 

^(ai + 29a2)-— ^ 



Cg = -(ai + 29a2)-— — (/3i+2«-i/J2), q = 2,3 



appearing in (47) are all equal to zero, given that we wish to maximise the order of the method. 
Thus, 

Co = a + ai + 1 = , 

Ci = ai + 2-(/3o + /3i+/32) =0, 

C2 = ^(ai+4)-(/3i+2/32)=0, 
Cs = ^(ai+8)-^(/3i+4/32)=0. 

Hence, 

ai = —1 — a , 

/3o = -^(l + 5a) , /3i = ^(l-a), /32 = ^(5 + a), 
and the resulting method is 

yn+2 - (1 + a)yn+i + ayn = [(5 + a)/„+2 + 8(1 - a)/„+i - (1 + 5a)/„] . (48) 

Further, 

C4 = ^(ai + 16)-l(/3i + 8/32) = -^(l + a) , 

C5 = |(ai + 32)-i(/?i + 16/32) = -^(17 + 13a). 

If a ^ — 1 then C4 7^ 0, and the method (48) is third order accurate. If, on the other hand, a = — 1, 
then C4 = and C5 ^ and the method (48) becomes Simpson's rule method - a fourth-order 
accurate two-step method. The error constant is: 

Ci = -^(1 + a), (49) 

C5 = a = -l. (50) 



Exercise 3 Determine all values of the real parameter b for which the linear multi-step 
method 

yn+3 + (26 - 3) (yn+2 - Vn+l) - Vn = hb{fn+2 + fn+l) 

is zero-stable. Show that there exists a value of b for which the order of the method is 4. 
Is the method convergent for this value ofb? Show further that if the method is zero-stable 
than its order cannot exceed 2. 
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Solution: According to the root condition, this hnear multi-step method is zero-stable if and 
only if all roots of its first characteristic polynomial 

p{z) = + i2b - 3){z^ - z) - I 

belong to the closed unit disc, and those on the unit circle arc simple. 

Clearly, p(l) = 0; upon dividing p{z) by z — 1 wc sec that p(z) can be written in the following 
factorised form: 

p{z) = {z-l) (z2 - 2(1 - b)z + l) = {z- l)pi(z) . 

Thus the method is zero-stable if and only if all roots of the polynomial pi {z) belong to the closed 
unit disc, and those on the unit circle are simple and differ from 1. Suppose that the method is 
zero-stable. Then, it follows that 6 7^ and 6 7^ 2, since these values of b correspond to double 
roots of pi{z) on the unit circle, respectively, z = 1 and z = —1. Since the product of the two roots 
of pi{z) is equal to 1 and neither of them is equal to ±1, it follows that they are strictly complex; 
hence the discriminant of the quadratic polynomial pi {z) is negative. Namely, 

4(1 - 6)2 - 4 < . 

In other words, b E (0, 2). 

Conversely, suppose that b G (0, 2). Then the roots of p{z) are 

01 = 1, Z2/3 = l-b + t^l-ib-lf . 

Since \z2/3\ = 1, 2:2/3 7^ 1 and 02 7^ 23, all roots of p{z) lie on the unit circle and they are simple. 
Hence the method is zero-stable. 

To summarise, the method is zero-stable if and only if 6 G (0, 2). 

In order to analyse the order of accuracy of the method we note that upon Taylor series 
expansion its truncation error can be written in the form 

Tn=(^l- ^) hY"{Xn) + J(6 - b)hV''{xn) + ^^(150 - 23b)h'y'' {Xn) + 0{h') . 

If 6 = 6, then T„ = 0(/i^) and so the method is of order 4. As 6 = 6 does not belong to the 
interval (0,2), we deduce that the method is not zero-stable for b = 6. 

Since zero-stability requires b G (0, 2), in which case 1 — g 7^ 0, it follows that if the method is 
zero-stable then Tn = 0{h'^). o 

3.4 Convergence 

The concepts of zero-stability and consistency are of great theoretical importance. How- 
ever, what matters most from the practical point of view is that the numerically computed 
approximations yn at the mesh-points x„, n = 0, . . . , N, are close to those of the analyt- 
ical solution y{xn) at these point, and that the global error e„ = y{xn) — yn between 
the numerical approximation ?/„ and the exact solution- value y{xn) decays when the step 
size h is reduced. In order to formalise the desired behaviour, we introduce the following 
definition. 

Definition 7 The linear multistep method (37) is said to be convergent if, for all initial 
value problems (1-2) subject to the hypotheses of Theorem 1, we have that 

lim yn = y{xn) (51) 

/i->0 

nh=x—xo 

holds for all x G [xo,-^m] md for all solutions {yn}n=o "/ difference equation (37) 
with consistent starting conditions, i.e. with starting conditions yg = r]s{h), s = 
0,1, . . . , k — 1, for which lim/j^o %(^) = Vo, s = 0,1, . . . , k — 1. 
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We emphasize here that Definition 7 requires that (51) hold not only for those sequences 
{yn}n=o which have been generated from (37) using exact starting vahics ys = y{xs), 
s = 0,1, . . . ,k — 1, but also for all sequences {yn}n=o whose starting values r]s{h) tend to 
the correct value, yo, as the h ^ 0. This assumption is made because in practice exact 
starting values are usually not available and have to be computed numerically. 

In the remainder of this section we shall investigate the interplay between zero-stability, 
consistency and convergence; the section culminates in Dahlquist's Equivalence Theorem 
which, under some technical assumptions, states that for a consistent linear multi-step 
method zero-stability is necessary and sufficient for convergence. 

3.4.1 Necessary conditions for convergence 

In this section we show that both zero-stability and consistency are necessary for conver- 
gence. 

Theorem 7 A necessary condition for the convergence of the linear multi-step method 
(37) is that it be zero-stable. 

Proof: Let us suppose that the linear multi-step method (37) is convergent; we wish to 
show that it is then zero-stable. 

We consider the initial value problem y' = 0, y(0) = 0, on the interval [0, Xm], Xm > 0, 
whose solution is, trivially, y{x) = 0. Applying (37) to this problem yields the difference 
equation 

OlkVn+k + Olk-lVn+k-l + ■ • • + OCoVn = . (52) 

Since the method is assumed to be convergent, for any a; > 0, we have that 

hm y„ = , (53) 

h^O 

nh=x 

for all solutions of (52) satisfying ys = ris{h), s = 0, . . . , fc — 1, where 

\mir]s{h)=Q, s = 0, 1, . . . A; - 1 . (54) 

Let z = re**^, be a root of the first characteristic polynomial p{z); r > 0, < < 2tt. It is 
an easy matter to verify then that the numbers 

yn = hr"' cos ncj) 

define a solution to (52) satisfying (54). If ^ 7^ and 4>i^i^, then 

Vl - yn+lVn-l _ ,2 2n 
sm (p 

Since the left-hand side of this identity converges to as /i — 0, n — 00, n/i = x, the 
same must be true of the right-hand side; therefore. 
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This implies that r < 1. In other words, we have proved that any root of the first 
characteristic polynomial of (37) lies in the closed unit disc. 

Next we prove that any root of the first characteristic polynomial of (37) that lies on 
the unit circle must be simple. Assume, instead, that z = re**^, is a multiple root of p{z); 
\r\ = 1, < 4> < 27r. We shall prove below that this contradicts our assumption that the 
method (52) is convergent. It is easy to check that the numbers 

Vn = h}^'^nr'^ cos{n(f)) (55) 

define a solution to (52) which satisfies (54), for 

\Vs{h)\ = \ys\ < h^'^s < h^l^(k - 1) , s = 0, . . . A; - 1 . 

If = or = TT, it follows from (55) with h = x/n that 

= . (56) 

Since, by assumption, l^l = 1, we deduce from (56) that lim|j__>Qo \yn\ = GO, which contra- 
dicts (53). 

If, on the other hand, / and / tt, then 

sin cj) 

where Zn = n^^h^^f^yn = h}f'^x~^yn- Since, by (53), lim„^oo -^n = 0, it follows that the 
left-hand side of (57) converges to as n ^ oo. But then the same must be true of the 
right-hand side of (57); however, the right-hand side of (57) of cannot converge to as 
n — >■ oo, since \z\ = 1. Thus, again, we have reached a contradiction. 

To summarise, we have proved that all roots of the first characteristic polynomial p{z) 
of the linear multi-step method (37) lie in the unit disc \z\ < 1, and those which belong 
to the unit circle = 1 are simple. By virtue of Theorem 6 the linear multi-step method 
is zero-stable. o 



Theorem 8 A necessary condition for the convergence of the linear multi-step method 
(37) is that it be consistent. 

Proof: Let us suppose that the linear multi-step method (37) is convergent; we wish to 
show that it is then consistent. 

Let us first show that Co = 0. We consider the initial value problem y' = 0, y{0) = 1, 
on the interval [0, Xm], Xm > 0, whose solution is, trivially, y{x) = 1. Applying (37) to 
this problem yields the difference equation 

otkVn+k + Oik-iUn+k-i + ■ • • + aoyn = . (58) 

We supply "exact" starting values for the numerical method; namely, we choose = 1, 
s = 0, . . . , A; — 1. Given that by hypothesis the method is convergent, we deduce that 

lim yn = l. (59) 

nh=x 
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Since in the present case y„ is independent of the choice of h, (59) is equivalent to saying 
that 

Hm yn = l. (60) 

n—>-oo 

Passing to the Umit n — oo in (58), we deduce that 

Ofe + ttk-i + . . . + ao = . (61) 

Recalling the definition of Co, (61) is equivalent to Co = (i.e. p(l) = 0). 

In order to show that Ci = 0, we now consider the initial value problem y' = 1, 
y(0) = 0, on the interval [0, Xm], Xm > 0, whose solution is y{x) = x. The difference 
equation (37) now becomes 

akVn+k + "fc-iyn+fc-i + • • ■ + aoVn = + Pk-l + • • ■ + /?o) , (62) 

where Xm — = Xm — = Nh and 1 < n < N — k. For a convergent method every 
solution of (62) satisfying 

lim?7s(/i) = , s = 0,l,...A;-l , (63) 
fe— >o 

where yg = r)s{h), s = 0,1, . . . ,k — 1, must also satisfy 

lim yn = X . (64) 

/i->0 

nh=x 

Since according to the previous theorem zero-stability is necessary for convergence, we 
may take it for granted that the first characteristic polynomial p{z) of the method does 
not have a multiple root on the unit circle \z\ = 1; therefore 

p'{l) = kak + ... + 2a2 + ai^0 . 

Let the sequence {yn}n=o be defined by y„ = Knh, where 

+ + + ; (65) 
kak + . . . + 2a2 + ai 

this sequence clearly satisfies (63) and is the solution of (62). Furthermore, (64) implies 
that 

X = yix) = lim ?/„ = lim Knh = Kx , 

/t->0 /t->0 

nh=x nh=x 

and therefore K = 1. Hence, from (65), 

Ci = {kak + . . . + 2a2 + ai) - (/3fe + . . . + /32 + /3i) = ; 
equivalently, p'(l) = cr(l). o 
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3.4.2 Sufficient conditions for convergence 

We begin by establishing some preliminary results. 

Lemma 1 Suppose that all roots of the polynomial p{z) = a^z^ + . . . + 012; + the 
closed unit disk 1^1 < 1 and those which lie on the unit circle \z\ = 1 are simple. Assume 
further that the numbers 7;, Z = 0, 1, 2, . . are defined by 

1 2 
rzr— r =70 + 712^ + 72^ +••• • 

Then, T = sup^>o 17(1 < 00 • 

Proof: Let us define p{z) = z^p{\/z) and note that, by virtue of our assumptions about 
the roots of p{z), the function 1/ p{z) is holomorphic in the open unit disc \z\ < 1. As the 
roots Zi, Z2, ■ ■ ■ , Zjn of p{z) on \z\ = 1 are simple, the same is true of the poles of l/p(z), 
and there exist constants Ai, . . ., A^ such that the function 

is holomorphic for \z\ < 1 and \f{z)\ < M for all \z\ < 1. Thus by Cauchy's estimate^ the 
coefficients of the Taylor expansion of / at z = also form a bounded sequence. As 

00 

Ai^zlz\ i = l,...,m, 
1=0 

and l^il < 1, we deduce from (66) that the coefficients in the Taylor series expansion of 
1/ p{z) form a bounded sequence, which completes the proof. o 

Now we shall apply Lemma 1 to estimate the solution of the linear difference equation 

Oik^m+k + «fe-iem+fe-l + • • ■ + OiQ^O = h{(5k,m.^m+k + Pk-l,m^m+k-l + ■ • • + Po,mem) + ■ 

(67) 

The result is stated in the next Lemma. 

Lemma 2 Suppose that all roots of the polynomial p{z) = a.kz^ + . . . + aiz + lie in the 
closed unit disk < 1 and those which lie on the unit circle \z\ = 1 are simple. Let B* 
and A denote nonnegative constants and (3 a positive constant such that 

\Pk,n\ + \h-l,n\ + ■■■ + |/3o,n| < -B* , 

|/3fe,nl</?' |An|<A, n = 0,l,...,Ar, 
and let < h < |ajfc|/3~^. Then every solution of (67) for which 

l^sl < E 1 s = 0, 1,...,A; — 1, 

^ Theorem (Cauchy's Estimate) If / is a holomorphic function in the open disc D(a, R), centre a 
and radius R, and \f{z)\ < M for all z G D{a, R), then \f^"\a)\ < M{n[/R"), n = 0, 1, 2, . . . . [For a proof 
of this result see, for example, Walter Rudin: Real and Complex Analysis. 3rd edition. McGraw-Hill, New 
York, 1986.] 
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satisfies 

\en\ < K* ex.p{nhL*) , n = (},l,...,N , 

where 

j^*^Y*B*, K* =r*{NA + AEk) , T* = T / {I - h\ak\~^ P) , 
r is as in Lemma 1, and 

A = Ittfcl + \ak-i\ + . . . + |ao| . 

Proof: For a fixed k we consider the numbers 7^, / = 0,1, . . . ,n — k, defined in Lemma 1. 
After multiplying both sides of the equation (67) for m = n — — Z by 7; , / = 0,1, ... ,n — k 
and summing the resulting equations, on denoting by Sn the sum, we find by manipulating 
the left-hand side in the sum that 

Sn = {oik^n + ak-ien-1 + . • • + Q;oe„_fe)7o 

+(Q!fce„-i + afe_ie„_2 + . . . + aQen-k-i)li + ■■■ 
+(afeefe + afe-iejt-i + . . . + Q;oeo)7n-fe 
= "fcToen + ("feTi + afe-i7o)en-i + • . . 

+(afe7n-fc + afe-i7n-fe-i + • • ■ + Q!o7n-2ifc)ejt 
+ {pLk-lln-k + • ■ • + Q;o7n-2fe+l)efe-l + . . . 
-|-ao7n-feeo • 

Defining 7; = for Z < and noting that 

f 1 for Z = 

(y-kll + otk-ili~i + ... + ao-fi-k = < Q for / > ^ ^ 

we have that 

Sn = en + {ak-ijn-k + • ■ • + ao7n-2fe+i)efe_i + . . . + ao7n-iteo ■ 
By manipulating similarly the right-hand side in the sum, we find that 

en + (afc_i7„_fc + . . . + ao7n-2fc+i)efc-i + • • • + ao7n-fceo 

= h [f3k,n-k70en + {I3k-l,n-k70 + f3k,n-k-l7l)en-l + ■■■ 
+ Wo,n-klO + • • • + f3k,n-2klk)en-k + • • • + Pofiln-k^o] 
+ (An-fc70 + A„_A:-l7l + • • • + Ao7n-fc) • 

Thus, by our assumptions and noting that by (68) 70 = a'^^, 

n—l 

Kl < hP\a^^\ \en\ + hVB* ^ |e^| + NTK + ATEk . 

m=0 

Hence, 

ji-i 

(1 - /i/3|a^^|)|en| < hVB* ^ |e„| + NTK + ATEk . 



m=0 
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Recalling the definitions of L* and K* we can rewrite the last inequality as follows: 

(t-i 

\en\<K* + hL*Y,\em\, n = 0, 1, . . . , TV . (69) 

m,=0 

The final estimate is deduced from (69) by induction. First, we note that by virtue of 
(68), AT>1. Consequently, K* > TAEk >Ek>E. Now, letting ra = 1 in (69), 

|ei| <K* + hL*\eo\ < K* + hL*E < K*{1 + hL*) . 

Repeating this argument, we find that 

\em\ < K^l + hL*)"^ , m = 0,...,fc-l. 

Now suppose that this inequality has already been shown to hold for m = 0, l,...,n— 1, 
where n > k; we shall prove that it then also holds for m = n, which will complete the 
induction. Indeed, we have from (69) that 

n -I- hr*\^ — 1 

Kl < K* + hL*K* ^^ I =K*{l + hL*Y . (70) 

Further, as 1 + hL* < e'*^* we have from (70) that 

|en| < if*e''^*" , n = 0,l,...,Ar . (71) 

That completes the proof of the lemma. We remark that the implication (69) (71) is 
usually referred to as the Discrete Gronwall Lemma. o 

Using Lemma 2 we can now show that zero-stability and consistency, which have been 
shown to be necessary are also sufficient conditions for convergence. 

Theorem 9 For a linear multi-step method that is consistent with the ordinary differen- 
tial equation (1) where f is assumed to satisfy a Lipschitz condition, and starting with 
consistent starting conditions, zero-stability is sufficient for convergence. 

Proof: Let us define 

6 = 6{h) = max |%(/i) — y{a + sh)\ ; 

0<s<fe— 1 

given that the starting conditions = r)s{h), s = 0, . . . ,k, are assumed to be consistent, 
we have that lim/i_>o S{h) = 0. We have to prove that 

lim y„ = y{x) 

n— »oo 

nh=x—xo 

for all X in the interval [xo,^m]- We begin the proof by estimating the truncation error 
of (37): 

(72) 



1 

ha{l) 



^ajy{xn+j) - h(3jy'{xn 



j=0 
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As y' G C[xo, Xm], it makes sense to define, for e > 0, the function 

X(e) = max \y'{x*) -y'{x)\ . 

X,X*e[xo,XM] 

For s = 0, 1, . . . , A;, we can then write 

y'{Xn+s) = y'{Xn) + 9sX{sh) , 

where l^^l < 1. Further, by the Mean Value Theorem, there exists £ {xn,Xn+s) such 
that 

y{Xn+s) = y{Xn) + Shy'{^s) ■ 

Thus, 

y{Xm+s) = y{Xm) + Sh [y' (Xm) + OgXish)] , 

where \dg\ < 1. 

Now we can write 

|o-(l)T„| < h~^{ai + a2 + . . . + ak)y{xn) + {ai + 2a2 + . . . + kak)y'{xn) 

- Wo + Pi + --- + Pk)y'{xn)\ 

+(|ai| + 2\a2\ + ... + k\ak\)\x{kh)\ + (|/3o| + \(3i\ + ... + m)\x{kh)\ . 

Since the method has been assumed consistent, the first, second and third terms on the 
right-had side cancel, giving 

|a(l)r„| < (|ai| + 2|q2| + . . . + k\ak\)\x{kh)\ + (|/?o| + + . . . + \(3k\)\x{kh)\ . 

Thus, 

\a{l)Tn\ < Kxikh) . (73) 

where 

K = \ai\ + 2\a2\ + ... + k\ak\ + \Po\ + + ... + \Pk\ ■ 
Comparing (37) with (72), we conclude that the global error e„i = y{xm) — ym satisfies 

"feem+ifc + ■ • ■ + "oeo = h {hdm+kf^m+k + • • • + f3o9mem) + cri^)Tnh , 

where 

^ _J [fiXm,y{Xm))-f{Xm,ym)]/e 

By virtue of (73), we then have that 

akCm+k + ■ • • + "oeo = h {hdm+kf^m+k + ■ • • + Pogmem) + dKx{kh)h . 

As / is assumed to satisfy the Lipschitz condition (3) we have that \gm\ < L, m = 0,1, . . .. 
On applying Lemma 2 with E = S{h), A = Kx{kh)h, N = {XM — xo)/h, B* = BL, where 
B = \l3o\ + l^il + . . . + \l3k\, we find that 

|e„| < r* [Ak6{h) + {Xm - xo)Kx{kh)\ exp[(x„ - xo)LT*B] , (74) 
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where 

r 

^ = l"o| + |ai| + • • • + |afe| , r* = rr^TTTTT ■ 

1 - h\af, f3k\L 

Now, y' is a continuous function on the closed interval [xq, Xm]', therefore it is uniformly 
continuous on [xo,Xm]- Thus, xi^h) — as ^ — 0; also, by virtue of the assumed 
consistency of the starting values, S{h) — as /t — 0. Passing to the limit /i — in (74), 
we deduce that 

lim le^l = ; 

n— *oo 

x—xo=nh 

equivalently, 

lim \y{x) - y„| = 

n—*oo 

x—xo=nh 

SO the method is convergent. o 

On combining Theorems 7, 8 and 9, we arrive at the following important result. 

Theorem 10 (Dahlquist) For a linear multi-step method that is consistent with the 
ordinary differential equation (1) where f is assumed to satisfy a Lipschitz condition, 
and starting with consistent initial data, zero-stability is necessary and sufficient for con- 
vergence. Moreover if the solution y{x) has continuous derivative of order {p + 1) and 
truncation error 0{hP), then the global error e„ = y{xn) — yn is also 0{hP). 

By virtue of Dahlquist's theorem, if a linear multi-step method is not zero-stable its 
global error cannot be made arbitrarily small by taking the mesh size h sufficiently small 
for any sufficiently accurate initial data. In fact, if the root condition is violated then 
there exists a solution to the linear multi-step method which will grow by an arbitrarily 
large factor in a fixed interval of x, however accurate the starting conditions are. This 
result highlights the importance of the concept of zero-stability and indicates its relevance 
in practical computations. 



3.5 Maximum order of a zero-stable linear multi-step method 

Let us suppose that we have already chosen the coefficients ctj, j = 0, . . . , /e, of the fe-step 
method (37). The question we shall be concerned with in this section is how to choose 
the coefficients Pj, j = 0, . . . ,k, so that the order of the resulting method (37) is as high 
as possible. 

In view of Theorem 10 we shall only be interested in consistent methods, so it is natural 
to assume that the ffist and second characteristic polynomials p{z) and a{z) associated 
with (37) satisfy p(l) = 0, p'{l) — a{l) = 0, with cr(l) / (the last condition being 
required for the sake of correctness of the definition of the truncation error (46)). 

Consider the function cp of the complex variable z, defined by 

Hz) = g - aiz) ; (75) 

the function logz appearing in the denominator is made single- valued by cutting the 
complex plane along the half-line ^z < 0. We begin our analysis with the following 
fundamental lemma. 
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Lemma 3 Suppose that p is a positive integer. The linear multistep method ( 37), with 
stability polynomials p{z) and cr(z), is of order of accuracy p if and only if the function 
4>{z) defined by (75) has a zero of multiplicity p at z = 1. 

Proof: Let us suppose that the fc-step method (37) for the numerical approximation of 
the initial value problem (1-2) is of order p. Then, for any sufficiently smooth function 
f{x,y), the resulting solution to (1-2) yields a truncation error of the form: 

r„ = %-/.V^+^)(x„) + o(/i^+i), 

(7(1) 

as /i — ^ 0, Cp+i ^ 0, Xn = xq + nh. In particular, for the initial value problem 

y' = y, yiO) = 1 , 

we get 



e 



nh 



p{e^') - ha{e^') = e"''%i/iP + 0{hP+^) , (76) 



ha{l) L"^' ' ' cr(l) 
as /t — 0, Cp+i 7^ 0. Thus, the function 

fih) = I [pie'') - haie'' 

is holomorphic in a neighbourhood of /i = and has a zero of order p at h = 0. The 
function z = e'^ is a bijective mapping of a neighbourhood oi h = onto a neighbourhood 
of 2; = 1. Therefore <p{z) is holomorphic in a neighbourhood of z = 1 and has a zero of 
multiplicity p at z = 1. 

Conversely, suppose that (j){z) has a zero of multiplicity p at z = 1. Then f{h) = (f>{e^) 
is a holomorphic function in the vicinity of /i = and has a zero of multiplicity path = Q. 
Therefore, 

k 

gih) = J2i<^,-hp,)e>'' 
3=0 

has a zero of multiplicity {p+ 1) at h = 0, implying that ^(0) = ^'(0) = . . . = ^^^^(0) = 0, 
but g''P+'^\0) / 0. First, 

k 

g{0) = = J2aj = Co. 

j=0 

Now, by successive differentiation of g with respect to h, 

k 

5'(0) = = ^(ia,-/3,) = Ci, 

j=0 
k 

g"{0)=0 = ^ifaj-2jP,)=2C2, 
j=0 

k 

g"'{0) = = J2if»J - 3//5i) = 2C3 , 
3=0 

etc. 

k 

g^Ho) = = ^(/a,- - pf-'Pj) =p\Cp. 
3=0 
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We deduce that Co = Ci = C2 = . . . = Cp = 0; since g^+^HO) 7^ we have that Cp+i 7^ 0. 
Consequently (37) is of order of accuracy p. o 

We shall use this lemma in the next theorem to supply a lower bound for the maximum 
order of a linear multistep method with prescribed first stability polynomial p{z). 

Theorem 11 Suppose that p(z) is a polynomial of degree k such that p{l) = and 7^ 
0, and let k he an integer such that < k < k. Then there exists a unique polynomial 
a{z) of degree k such that p'{l) — cr(l) = and the order of the linear multi-step method 
associated with p{z) and a{z) is > k + 1. 

Proof: Since the function p{z)/log{z) is holomorphic in the neighbourhood of z = 1 it 
can be expanded into a convergent Taylor series: 

f^=C0+Cl{z-l) + C2{z-lf + ... . 

logz 

On multiplying both sides by log^; and differentiating we deduce that cq = (7^ 0). 
Let us define 

a{z) = Co + ci(z - 1) + . . . + Cj^iz - 1)^= . 
Clearly a{l) = cq = p'(l) (7^ 0). With this definition, 



and therefore ^(z) has a zero at 2; = 1 of multiplicity not less than k + 1. By Lemma 3 
the linear fc-step method associated with p{z) and a{z) is of order > k + 1. 

The uniqueness of (t(z) possessing the desired properties follows from the uniqueness 
of the Taylor series expansion of (p(z) about the point z = 1. o 

We note in connection with this theorem that for most methods of practical interest 
either k = k — 1 resulting in an explicit method or k = k corresponding to an implicit 
method. In the next example we shall encounter the latter situation. 

Example 5 Consider a linear two-step method with p{z) = {z — \){z — X) . The method 
will he zero-stable as long as X £ [—1,1). Consider the Taylor series expansion of the 
function p{z)/log{z) ahout the point z = 1: 

pjz) ^ jz-m-x+jz-i)) 

logz log[l + (z- 1)] 

- [l-A + (.-l)lx{l-^+<^-<^ + 0((.-lrtp 

= i-A + ^(.-i) + ^(.-i)2-i±^(.-iy'+o((.-i)'). 

A two-step method of maximum order is obtained by selecting 

a{z) = l-A + ^(.-l) + ^(.-l)^ 
1 + 5A 2-2A 5 + A 2 
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If —1, the resulting method is of third order, with error constant 



1 + A 

^^ = "^' 

whereas if X = —1 the method is of fourth order. 
In the former case the method is 

yn+2 -{l + X) + Xyn = h \—^fn+2 + ^ fn+\ ) 

with X a parameter contained in the interval (—1,1). In the latter case, the method has 
the form 

yn+2 -yn= 2 (/n+2 + 4/„+i + /„) , 

and is referred to as Simpson's method. 

By inspection, the linear fc-step method (37) has 2k+2 coefficients: aj, Pj, j = 0, . . . , k, 
of which ak is taken to be 1 by normaUsation. This leaves us with 2k + 1 free parameters 
if the method is implicit and 2k free parameters if the method is implicit (given that in 
the latter case (dk is fixed to have value 0). According to (47), if the method is required to 
have order p, p + 1 linear relationships Cq = 0, . . . ,Cp = involving aj, (dj, j = 0, . . . , k, 
must be satisfied. Thus, in the case of the implicit method, we can impose p + 1 = 2k + l 
linear constraints Co = 0, . . . , C2k+i = to determine the unknown constants, yielding 
a method of order p = 2k. Similarly, in the case of an explicit method, the highest order 
we can expect is p = 2k — I. Unfortunately, there is no guarantee that such methods will 
be zero-stable. Indeed, in a paper published in 1956 Dahlquist proved that there is no 
consistent, zero-stable A;-step method which is of order > {k + 2). Therefore the maximum 
orders 2k and 2k — 1 cannot be attained without violating the condition of zero-stability. 
We formalise these facts in the next theorem. 

Theorem 12 There is no zero-stable linear k-step method whose order exceeds k + 1 if k 
is odd or k + 2 if k is even. 

Proof: Consider a linear A;-step method (37) with associated first and second stability 
polynomials p and a. Further, consider the transformation 

C + 1 

which maps the open unit disc |C| < 1 of the ^-plane onto the open half-plane < of 
the z-plane; the circle |C| = 1 is mapped onto the imaginary axis = 0, the point C = 1 
onto z = 0, and the point C = — 1 onto ( = oo. 
It is clear that the functions r and s defined by 



.2 7'^ VI -2/ \ 2 J \1- z 

are in fact polynomials, deg(r) < k and deg(s) < k. 
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If p{Q has a root of multiplicity p, < p < k, at ( = (^q ^ —1, then r{z) has a root 
of the same multiplicity at z = (Co — l)/(Co + 1); if p(C) ^ root of multiplicity p > I, 
< p < k, at ( = —1, then r{z) is of degree k — p. 

Since, by assumption, the method is zero-stable, C = 1 is a simple root of p(C); conse- 
quently, z = is a simple root of r{z). Thus, 

r{z) = aiz + a2Z^ + . . . -|- akZ^ , ai / , G R . 

It can be assumed, without loss of generality, that ai > 0. Since by zero stability all roots 
of p{() are contained in the closed unit disc, we deduce that all roots of r{z) have real 
parts < 0. Therefore, all coefficients aj, j = 1, . . . ,k, of r{z) are nonnegative. 
Now let us consider the function 

i-z\'' , ri + z\ 1 



Qiz) = [-^) <^ j = r{z) - s{z) . 

The function q{z) has a zero of multiplicity p at 2; = if and only if defined by (75) 

has a zero of multiplicity p at ( = I; according to Lemma 3 this is equivalent to the linear 
k-step method associated with p{() and a{0 having order p. Thus if the linear fc-step 
method associated with p{z) and a{z) has order p then 

s{z) = bo + biz + b2Z^ + ... + bp-izP~^ , 

where 

1x7 — = bo + biz + b2Z^ + ... . 



As the degree of s{z) is < k, the existence of a consistent zero-stable k-step linear multistep 
method of order p> k + 1 {ot p> k + 2) now hinges on the possibility that 

bk+i = . . . = bp-i = , (or bk+2 = ■■■ = = 0) . 

Let us consider whether this is possible. 

We denote by cq, ci, C2, . . ., the coefficients in the Taylor series expansion of the function 



logi 



l+z 

z 



namely. 



-j^ = Co + C2Z'^ + C4Z^ + ... . 



Then, adopting the notational convention that a^, = loi v > k, we have that 

60 = coao , 

bi = coa2 , 
etc. 

b2u = coa2i/+i + C2a2u-i + • • • + C2^ai , 

b2u+l = Co02z/+2 + C2a2i/ + . . . + C2i/a2 , z/=l,2, .... 

It is a straightforward matter to prove that C2u < 0, v = 1,2, . . . (see also Lemma 5.4 on 
page p. 233 of Henrici's book). 
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(i) If k is an odd number, then, since 0,^ = for > A;, we have that 

h+i = C2ak + C4afc-2 + • • • + Cfc+iai . 
Since ai > and no ai, is negative, it follows that b^+i < 0. 

(ii) If A; is an even number, then 

^fc+i = C2ak + Ciak-2 + • • • + Cfca2 . 

Since C2v < 0, = 1, 2, . . . , and > 0, fi = 2,3, . . . , k, we deduce that b^+i = if 
and only if 02 = 04 = . . . = = 0, i.e. when r(z) is an odd function of z. This, 
together with the fact that all roots of r{z) have real part < 0, implies that all roots 
of r(z) mush have real part equal to zero. Consequently, all roots of p{C) lie on 
Id = 1. Since = 0, the degree of r{z) is A; — 1, and therefore —1 is a (simple) root 
ofp(C). 

As C2v < 0, > and ai > 0, it follows that 

bk+2 = Cittk-i + ceUk-s + • • ■ + Ck+2ai < , 

showing that bk+2 7^ 0- 

Thus if a fc-stcp method is zero-stable and k is odd then bk+i 7^ 0, whereas if k is even 
then bk+2- This proves that there is no zero-stable fc-step method whose order exceeds 
A + 1 if A; is odd or A; + 2 if A is even. o 

Definition 8 A zero-stable linear k-step method of order k + 2 is said to be an optimal 
method. 

According to the proof of the previous theorem, all roots of the first characteristic 
polynomial p associated with an optimal linear multistep method have modulus 1. 

Example 6 As k + 2 = 2k if an only ifk = 2 and Simpson's rule method is the zero-stable 
linear 2-step method of maximum order, we deduce that Simpson's rule method is the only 
zero-stable linear multistep method which is both of maximum order {2k = 4) and optimal 
(A + 2 = 4). 

Optimal methods have certain disadvantages in terms of their stability properties; we 
shall return to this question later on in the notes. 

Linear A-step methods for which the first characteristic polynomial has the form 
p{z) = z^ — z^~^ are called Adams methods. Explicit Adams methods are referred to 
as Adams— Bashforth methods, while implicit Adams methods are termed Adams— 
Moulton methods. Linear A-step methods for which p{z) = z^ — z^~'^ are called 
Nystrom methods if explicit and Milne— Simpson methods if implicit. All these 
methods are zero-stable. 
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3.6 Absolute stability of linear multistep methods 

Up to now we have been concerned with the stability and accuracy properties of linear 
multistep methods in the asymptotic limit of /i ^ 0, n ^ oo, nh fixed. However, it is of 
practical significance to investigate the performance of methods in the case of /i > fixed 
and n ^ oo. Specifically, we would like to ensure that when applied to an initial value 
problem whose solution decays to zero as a; — oo, the linear multistep method exhibits a 
similar behaviour, for h > fixed and Xn = Xq + nh ^ oo. 

The canonical model problem with exponentially decaying solution is 

y' = Xy, x>0, y(0) = yo (t^ 0) , (77) 

where 3f?A < 0. Indeed, 

y{x) = yoe'^^^^e"^^ , 

and therefore, 

\y{x)\ < cxp(-.T|KA|) , .T > , 

yielding lim^^^oo y{x) = . Thus, using the terminology introduced in the last paragraph 
of Section 1, the solution is asymptotically stable. 

In the rest of the section we shall assume, for simplicity, that A is a negative real 
number, but everything we shall say extends straightforwardly to the general case of A 
complex, 5RA < 0. 

Now consider the linear fc-step method (37) and apply it to the model problem (77) 
with A real and negative. This yields the linear difference equation 

k 

{aj - h\(3j)yn+j = . 

j=0 

Given that the general sohition y„ to this homogeneous difference equation can be ex- 
pressed as a linear combination of powers of roots of the associated characteristic polyno- 
mial 

tt{z; h) = p{z) - ha{z) , {h = h\) , (78) 

it follows that y„ will converge to zero for h > fixed and — >■ oo if and only if all roots of 
tt{z; h) have modulus < 1. The A;th degree polynomial ^{z] h) defined by (78) is called the 
stability polynomial of the linear /c-step method with first and second characteristic 
polynomials p{z) and cj(z), respectively. This motivates the following definition. 

Definition 9 The linear multistep method (37) is called absolutely stable for a given 
h if and only if for that h all the roots rg = rsih) of the stability polynomial 7r(z, h) defined 
by (78) satisfy \rs\ < I, s = l,...,k. Otherwise, the method is said to be absolutely 
unstable. An interval (a, f3) of the real line is called the interval of absolute stability 
if the method is absolutely stable for all h € {a, 13). If the method is absolutely unstable 
for all h, it is said to have no interval of absolute stability. 

Since for A > the solution of (77) exhibits exponential growth, it is reasonable 
to expect that a consistent and zero-stable (and, therefore, convergent) linear multistep 
method will have a similar behaviour for /t > sufficiently small, and will be therefore 
absolutely unstable for small h = Xh. According to the next theorem, this is indeed the 
case. 
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Theorem 13 Every consistent and zero-stable linear multistep method is absolutely un- 
stable for small positive h. 

Proof: Given that the method is consistent, there exists an integer p > 1 such that 
Co = Ci = . . . = Cp = and Cp+i ^ 0. Let us consider 



7r(e^;/i) = p(e^) - /i(7(e") 

k 



j=0 
k 

E 

k 

E 

j=0 

k 



oo 



h(3je' 



hj 



9=0 

oo 



q=0 

oo 



«. E , 

9=1 ^ 

k 



9=0 

= E"i + E 

j=0 j=0 

j=0 9=1 

oo 

9=1 

oo 
9=P+1 



E" 

j=0 



%^i(^-i)' 

fc ^-9-1 



(79) 



On the other hand, noting that the polynomial Tr{z; h) can be written in the factorised 
form 

it{z, h) = {ak - hPk){z -ri)...{z- r^) 
where = rgih), s = 1, . . . , k, signify the roots of 7r(.; h), we deduce that 

7r(e^ h) = {ak - hf3k){e'' - ri{h)) ...{e^- rk{h)) . (80) 

As ^ — 0, ajk — h/3k — > ctfe 7^ and rs(h) s = 1, . . . ,k, where s = 1, . . . ,k, 

arc the roots of the first stability polynomial C{^)- Since, by assumption, the method is 
consistent, 1 is a root of C{^)i furthermore, by zero-stability 1 is a simple root of C{^)- Let 
us suppose, for the sake of definiteness that it is (i that is equal to 1. Then, Cs 7^ 1 for 
s 7^ 1 and therefore 

lim(e'^-r,(/i)) = (l-Cs)7^0, 

We conclude from (80) that the only factor of Tr{e'^;h) that converges to as ^ is 
— ri(h) (the other factors converge to nonzero constants). Now, by (79), Tr{e^;h) = 



0{hP+^), so it follows that 



- ri{h) = 0{hP+^) . 
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Thus we have shown that 

ri{h) = + 0(hP+^) . 

This imphes that 

ri(h) > 1 + 

for small positive h. That completes the proof. o 

According to the definition adopted in the previous section, an optimal A;-step method 
is a zero-stable linear /c-step method of order k + 2. We have also seen in the proof 
of Theorem 12 that all roots of the first characteristic polynomial of an optimal /c-step 
method lie on the unit circle. By refining the proof of Theorem 13 it can be shown that 
an optimal linear multistep method has no interval of absolute stability. 

It also follows from Theorem 13 that whenever a consistent zero-stable linear multistep 
method is used for the numerical solution of the initial value problem (1-2) where > 0, 
the error of the method will increase as the computation proceeds. 

3.7 General methods for locating the interval of absolute stability 

In this section we shall describe two methods for identifying the endpoints of the interval 
of absolute stability. The first of these is based on the Schur criterion, the second on the 
Routh-Hurwitz criterion. 

3.7.1 The Schur criterion 

Consider the polynomial 

(j){r) = Ckv'^ + . . . + cir + Co , ^ , cq 7^ , 

with complex coefficients. The polynomial (p is said to be a Schur polynomial if each 
of its roots Vg satisfies jr^l < 1, s = 1, . . . , k. 
Let us consider the polynomial 

(t){r) = Cor'' + cir^-^ + . . . + c^-ir + Ck , 

where Cj denotes the complex conjugate of Cj, j = 1, . . . ,k. Further, let us define 

'm<Pir)-mHr) ■ 

Clearly has degree < A; — 1. 

The following key result is stated without proof. 

Theorem 14 (Schur's Criterion) The polynomial (j) is a Schur polynomial if and only 
i/ 10(0)1 > |<?!>(0)| and 4>i is a Schur polynomial. 

We illustrate Schur's criterion by a simple example. 

Exercise 4 Use Schur's criterion to determine the interval of absolute stability of the 
linear multistep method 

yn+2 -yn = (/n+1 + 3/n) . 



Mr) = - 
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Solution: The first and second characteristic polynomials of the method are 

1 



p{z) =z'-l, 
Therefore the stability polynomial is 



^{z) = -(z + 3) 



7r(r; h) = p{r) — ha{r) = 



-hr ■ 



Now, 



7t{r; h) 



3-\ 1- 

-h - -hr + l . 

2 2 



Clearly, |7r(0;/i)| > |7r(0,/i)| if and only if /i e (-|,0). As 

has the imiquc root — | and is, therefore, a Schur polynomial, we deduce from Schur's criterion 
that Tr{r;h) is a Schur polynomial if and only ii h G (— 1,0). Therefore the interval of absolute 
stability is (— 1,0). o 

3.7.2 The Routh— Hurwitz criterion 

Consider the mapping 



1 



z = 



r+1 

of the open unit disc |r| < 1 of the complex r-plane to the open left half-plane < of 
the complex z-plane. The inverse of this mapping is given by 

l + z 

r = . 

1 — z 

Under this transformation the function 



becomes 



7r(r, h) = p{r) — ha{r) 



'l + z\ - fl + z 
p 1 z I - ha 



l-zj \l-z. 
On multiplying this by (1 — z)'^, we obtain a polynomial of the form 

aoz'' + aiz^~'^ + . . . + Ofe . 



(81) 



The roots of the stability polynomial 7r(r, h) lie inside the open unit disk |r| < 1 if and 
only if the roots of the polynomial (81) lie in the open left half-plane 3f?z < 0. 

Theorem 15 (Routh— Hurwitz Criterion) The roots of (81) lie in the open left half- 
plane if and only if all the leading principal minors of the k x k matrix 



Q 



ai 


as 


05 . 


• • a2fc~i 


0,0 


0-2 


O4 


• • a2fc-2 





Oi 


03 . 


• • a2fc-3 





ao 


02 . 


• • a2fe-4 








. 





are positive and ao > 0; we assume that aj = if j > k. In particular, 



46 



a) for k = 2: oq > 0, ai > 0, 02 > 0, 

b) for k = 3; ao > 0, oi > 0, 02 > 0, 03 > 0, 0102 — 0300 > 0, 

c) for k = 4; oq > 0, ai > 0, 02 > 0, > 0, 04 > 0, 010203 — 00O3 — 040^ > 

represent the necessary and sufficient conditions for ensuring that all roots of (81) lie in 
the open left half-plane. 

We illustrate this result by a simple exercise. 

Exercise 5 Use the Routh-Hurwitz criterion to determine the interval of absolute stability 
of the linear multistep method from the previous exercise. 

Solution: On applying the substitution 

l + z 
r = 



in the stability polynomial 

7r(r, h) = r^ — ^hr ~ 
and multiplying the resulting function by (1 — 2)^, we get 



(1-^) 

with 



i-z 2 vi-^y V 2 



= aoz^ + aiz + 02 



flo = —h , fli = 4 + 3/1 , 0,2 = —2h 



Applying part a) of Theorem 15 we deduce that the method is zero-stable if and only ii h £ (— |, 0); 
hence the interval of absolute stability is (— 1,0). o 

Wc conclude this section by listing the intervals of absolute stability (a, 0) of /c-step 
Adams-Bashforth and Adams-Moulton methods, for k = 1,2,3,4. We shall also supply 
the orders p* and p and error constants Cp*+i and Cp+i, respectively, of these methods. 
The verification of the stated properties is left to the reader as exercise. 

A;-step Adams— Bashforth (explicit) methods: 

(1) k = l,p* = l,Cp*+i = l,a = -2, 

yi-yo = hfo ; 

(2) k = 2,p* = 2, Cp*+i = ^,a = -l, 

y2-yi = ^(3/1 - /o) ; 

(3) k = 3,p* = 3,Cp*+i = la = -^, 



ys- 2/2 = ^(23/2-16/1 + 5/0); 
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(4) k = 4,p* = A, Cp*+i = is, a = , 



y4 - 2/3 = ^(55/3 - 59/2 + 37/i - 9/o) . 



fc-step Adams— Moulton (implicit) methods: 

(1) k = l,p = 2, Cp+i = -j^, a = -oo , 

yi-yo = ^(/i + /o) ; 



(2) k = 2,p = 3, Cp+1 = -^,a = -6, 



y2-yi = ^i^h + 8/i - /o) ; 



(3) k = 3,p = 4, Cp+i = a = -3 , 



2/3 - y2 = ^(9/3 + 19/2 -5/1 + /0); 



(4) k = 4,p = 5, Cp+i = 



1440' " 49 ' 



2/4 - 2/3 = ^(251/4 + 646/3 - 264/2 + IO6/1 - 19/o) . 

We notice that the A;-step Adams-Moulton (implicit) method has larger interval of absolute 
stability and smaller error constant than the fc-step Adams-Bashforth (explicit) method. 

3.8 Predictor-corrector methods 

Let us suppose that we wish to use the implicit linear fc-step method 

k k 

X] ^jVn+j = ^ X] Pjfn+j , Oik, (3k ■ 

j=0 j=0 

Then, at each step we have to solve for yn+k the equation 

fc-i 

Olkyn+k - h(3kf{Xn+k,yn+k) = X (^(^jfn+j " Oljyn+j) ■ 

3=0 

li h < \ak\l {L\l3k\) where L is the Lipschitz constant of / with respect to y (as in Pi- 
card's Theorem 1), then this equation has a unique solution, yn+k'-, moreover, yn+k can be 
computed by means of the fixed-point iteration 

k-l k-l 

(Xkyntk + Yl "jVn+j = hPkf{Xn+k, yi+fe) + ^ X ' S = 1, 2, 3, . . . , 

j+0 j=0 
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with y^^]_j^ a suitably chosen starting value. 

Theoretically, we would iterate until the iterates have converged (in practice, until 

some stopping criterion such as — y^+kl < ^ is satisfied, where e is some preassigned 

tolerance). We would then regard the converged value as an acceptable approximation 
Un+k to the unknown analytical solution- value y{xn+k)- This procedure is usually referred 
to as correcting to convergence. 

Unfortunately, in practice, such an approach is usually unacceptable due to the amount 
of work involved: each step of the iteration involves an evaluation of /(x^+fe, y|*^^,) which 
may be quite time-consuming. In order to keep to a minimum the number of times 
f{xn+k^yn\-k) evaluated, the initial guess y^+k must be chosen accurately. This is 
achieved by evaluating yll^j.^ by a separate explicit method called the predictor, and 
taking this as the initial guess for the iteration based on the implicit method. The implicit 
method is called the corrector. 

For the sake of simplicity we shall suppose that the predictor and the corrector have the 
same number of steps, say k, but in the case of the corrector we shall allow that both ao 
and Po vanish. Let the linear multistep method used as predictor have the characteristic 
polynomials 

p*{z) = J2a*z^ , 4 = 1, a*{z) = ^P*z^, (82) 

i=o j=Q 

and suppose that the corrector has characteristic polynomials 

k k 

P{^) = J2 "i^^ ' "fc = 1 , (^{z) = ^ PjZ^ . (83) 

3=0 j=0 

Suppose that m is a positive integer: it will denote the number of times the corrector is 
applied; in practice m < 2. If P indicates the application of the predictor, C a single 
application of the corrector, and E an evaluation of / in terms of the known values of its 
arguments, then P{EC)"^E and P{EC)"^ denote the following procedures. 

a) P{EC)'^E 

fc-1 fc-1 

Vnlk + ' 
j=0 



fe-1 



j=0 



yn+j 






j=0 


As] 
J n+k 


= f{Xn+k, ylXk) ' 






H 






i=o 


Am] 
J n+k 


= f{Xn+ki yi+fe) ) 



s = 0, . . . , m — 1, 



forn = 0,1,2,.... 

b) PiEC)"" 

fe-1 fc-1 

Vn+k + "jWn-hj - Z^l^j Jn+j ' 
j=0 j=0 
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fc-1 

iln+k ^ 

j=0 

for n = 0,1,2,.. 



Jn+k 



[m\ 



fc-1 



3=0 



0, . . . , m — 1 , 



3.8.1 Absolute stability of predictor-corrector methods 

Let us apply the predictor-corrector method P{EC)"^E to the model problem 

y' = Xy, y{0) = yo 0) , 



(84) 



where A < 0, whose solution is, trivially, the decaying exponential y{x) = yocxp{Xx), 
X > 0. Our aim is to identify the values of the step size h for which the numerical solution 
computed by the P{EC)'^E method exhibits a similar exponential decay. The resulting 
recursion is 



fc-1 



fc-1 



i/n+fc Vn+j — Z^Pj Vn+j ' 

y?S' + E«.-yH, = hpJX, + -hY.P,y^-},^ s = 0,...,m-l, 

j=o j=0 

for n = 0, 1,2, . . ., where h = \h. In order to rewrite this set of equations as a single 
difference equation involving yll"^ , , . . . y^^+k O'^tyj have to eliminate the intermediate 
stages involving • • • , y^n+k^ from the above recursion. 

Let us first take s = and eliminate y^^j^}. form the resulting pair of equations to obtain 



fc-1 



fc-1 



fc-1 



fe-i 



^Sfc + E = h E - E + ^ E P^n. 



n+j 



3=0 



j=0 



j=0 



3=0 



Now take s = 1 and use the last equation to eliminate y|j^^j.; this gives, 



fc-1 

E' 

3=0 



[2] \ - [mj 
Vn+k ^jyn+j 



hpk 



fc-1 



fc-1 



fc-1 



fc-1 



j=0 j=0 



fc-1 



+ ^E/5iyH. 



Equivalently, 



^ifc + (i + WE«.ya 

i=o 

= f;^E - E + (1 + ^/?^)^ E /5.*'j • 

V :/=o 3=0 J j=o 
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By induction, 

fe-i 



3=0 

I fe— 1 k—\ \ k—\ 



mj 

j=0 j=0 ) ' j=0 

For m fixed, this is a kth. order difference equation involving . . . , yj^^- In order to 
ensure that the solution to this exhibits exponential decay as n — ^ oo, we have to assume 
that all roots to the associated characteristic equation 

fc-i 



j=o 

k-l k-1 \ fc-1 

X ^ ^ J; l' X ^ .1-1.1/^ . / 1 \ ■il l 1 \ 1 \. ^ ,1 

j=0 j=0 I j=0 



have modulus < 1. This can be rewritten in the equivalent form 

Az'' + {i+h(3k+... + mr-') ip{z) - ha{z)) + mr - h^*i^)) = o , 

where 

A = i + (i+hPk + ...+ mr~') m - + mrm - «d , 

Now, since ccfe = = 1 and = 0, we deduce that ^ = 0, and therefore the characteristic 
equation of the P{ECrE method is 

t^p{ec)^e{z; h) = p{z) - ha{z) + M„^{h){p*{z) - ha*{z)) = , 

where 

""W = i + ftA + ... + (Aa)..-i . "-SI. 

Here, T^p[EC)^E{^'ih) is referred to as the stability polynomial of the predictor-corrector 
method P{ECrE. 

By pursuing a similar argument we can also deduce that the characteristic equation of 
the predictor corrector method P[ECr is 

■^p{EC)r-{z\ h) = p{z) - hu{z) + ^^^^^ {P*{z)(^{z) - p{z)a*{z)) = . 

Here, TTp(_EC)"^iz',h) is referred to as the stability polynomial of the predictor-corrector 
method P{ECr- 

With the predictor and corrector specified, one can now check using the Schur criterion 
or the Routh-Hurwitz criterion, just as in the case of a single multi-step method, whether 
the roots of 'Kp(^^c^m^{z] h) and 7rp(g(^-)m (z; h) all lie in the open unit disk \z\ < 1 thereby 
ensuring the absolute stability of the P{ECrE and P{ECr method, respectively. 
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Let us suppose, for example, that \hPk\ < 1, i.e. that < ^ < l/|A/3jfc|; then, 
hm^_»oo Mm(h) = 0, and consequently, 

lim ttp(ec)^e{z\ h) = tt{z, h) , lim 7rp(EC)-(^; h) = it{z, h) , 

where 7r(z; h) = p{z) — ha{z) is the stability polynomial of the corrector. This implies 
that in the mode of correcting to convergence the absolute stability properties of the 
predictor-corrector method are those of the corrector alone, provided that \hPk\ < 1. 



3.8.2 The accuracy of predictor-corrector methods 

Let us suppose that the predictor P has order of accuracy p* and the corrector has order of 
accuracy p. The question we would like to investigate here is: What is the overall accuracy 
of the predictor- corrector method? 

Let us consider the P{EC)^E method applied to the model problem (84) with m > 1. 
We have that 

7rp(EC)-E(c^ = />(c^) - M^^) + MM{p*{e^) - ha*{e^)) 
= 0{K^+') + M^ih)OihP*+^) 
= 0{hP+'^ + 

0(hP+^ + hP+^) iip*>p 

0(hP~^^) if p* = p — q, < q < p and m> q 

0{hP+^ + hP-i+"'+^) ifp* =p-q,0< q<p and m<q-l. 

Consequently, denoting by T„ ^ ' the truncation error of the method P{EC)"'-E, we 
have that 



rpP{EC)'^E 



0{hP) iip*>p 

0{hP) ii p* = p — q, < q < p and m > q 

0{hP-i+'^) ifp* = p-q,0<q<pandm<q-l. 



This implies that from the point of view of overall accuracy there is no particular advantage 
in using a predictor of order p* > p. Indeed, as long as p* + m > p, the predictor-corrector 
method P(EC)^E will have order of accuracy p. 

Similar statements can be made about P{EC)"^ type predictor-corrector methods with 
m > 1. On writing 

p*{z)a{z) - a*{z)p{z) = {p*{z) - ha\z))a{z) - a*{z){p{z) - ha{z)) , 

we deduce that 

7rp(EC)-(c^ h) = 0{hP+^ + hP'+"' + hP+"') . 

Consequently, as long as p* -|-m >p-\-l the predictor-corrector method P{EC)"^ will have 
order of accuracy p. 
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4 Stiff problems 



Let us consider an initial value problem for a system of m ordinary differential equations 
of the form: 

y' = f{x,y), y(a) = yo, (85) 

where y = (yi, . . . ,ym)'^- A linear A;-step method for the numerical solution of (85) has 
the form 

k k 

^jyn+j = hY, ^j^ri+j , (86) 

3=0 3=0 

where in+j = ^{'^n+jiVn+j)- Let us suppose, for simplicity, that f(x,y) = ^y + b where 
^ is a constant matrix of size mx m and b is a constant (column) vector of size m; then 
(86) becomes 

k 

^(q,/ - hp,A)yn+j = ha{l)h , (87) 

j=0 

where cr{l) = J2j=o Pj (7^ 0) and / is the m x rn identity matrix. Let us suppose that the 
eigenvalues Aj, i = 1, . . . , m, of the matrix A are distinct. Then, there exists a nonsingular 
matrix H such that 



HAH~^ = A = 



Ai 


. 


.. 





A2 . 


.. 





. 


■ • -^m 



Let us define z = Hy and c = ha{l)Hh. Thus, (87) can be rewritten as 

k 

Yiajl - hPjA)zn+j = c , (88) 
3=0 

or, in component-wise form, 

k 

^(aj - hf3j\i)zn+j,i = Ci , 
3=0 

where Zn+j,i and Cj, i = 1, . . . ,m, are the components of z^+j and c respectively. Each 
of these m equations is completely decoupled from the other m — 1 equations. Thus we 
are now in the framework of Section 3 where we considered linear multistep methods for a 
single differential equation. However, there is a new feature here: given that the numbers 
Aj, i = 1, . . . ,m, are eigenvalues of the matrix A, they need not be real numbers. As a 
consequence the parameter h = hX, where A is any of the m eigenvalues, can be a complex 
number. This leads to the following modification of the definition of absolute stability. 

Definition 10 A linear k-step method is said to be absolutely stable in an open set TZa 
of the complex plane, if for all h G TZa all roots rg, s = 1, . . . , k, of the stability polynomial 
Tr{r,h) associated with the method, and defined by (78), satisfy \rs\ < 1. The set TZa is 
called the region of absolute stability of the method. 

Clearly, the interval of absolute stability of a linear multistep method is a subset of its 
region of absolute stability. 
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Exercise 6 

a) Find the region of absolute stability of Euler's explicit method when applied to the 
ordinary differential equation y' = Xy, y{xo) = yg. 

b) Suppose that Euler's explicit method is applied to the second-order differential equa- 
tion 

y" + (A + l)y' + Ay = , y{0) = 1 , y'(0) = A - 2 , 

rewritten as a first-order system in the vector [u, v)'^ , with u = y and v = y' . Suppose 
that A >> 1. What choice of the step size h will guarantee absolute stability in the 
sense of Definition 10? 

Solution: a) For Euler's explicit method p{z) = z — 1 and a{z) = 1, so that ^{z; h) = p{z) — 
ha{z) = {z — 1) — h = z — {1 + h). This has the root r = (1 + h). Hence the region of absolute 
stability is TZa = {h & C : \1 -\-h\ < 1} which is an open unit circle centred at —1. 

b) Now writing u = y and v = y' and y = (u, w)-^, the initial value problem for the given 
second-order differential equation can be recast as 

y' = Ay, y(0) = yo , 

where 

^=(-A -(A + 1)) ^'^^ y'={x-2) ■ 
The eigenvalues of A are the roots of the characteristic polynomial of A, 

det{A - zl) = z^ + {\ + l)z-^\ . 

Hence, ri = —1, r2 = —A; we deduce that the method is absolutely stable provided that h <2/X. 
It is an easy matter to show that 

u{x) = 2e-^ - e"^^ , v{x) = -2e-^ + Ae"^^ . 

The graphs of the functions u and v arc shown in the figure below for A = 45. Note that v exhibits 
a rapidly varying behaviour (fast time scale) near x = while u is slowly varying for a; > and 

V is slowly varying for x > 1/45. Despite the fact that over most of the interval [0, 1] both u and 

V are slowly varying when A — 45, we are forced to use a step size of ft. < 2/45 in order to ensure 
that the method is absolutely stable. o 

In the example the v component of the solution exhibited two vastly different time 
scales; in addition, the fast time scale (which occurs between x = and x « 1/A) has 
negligible effect on the solution so its accurate resolution does not appear to be important 
for obtaining an overall accurate solution. Nevertheless, in order to ensure the stability 
of the numerical method under consideration, the mesh size h is forced to be exceedingly 
small, h < 2/X, smaller than an accurate approximation of the solution for x 3> 1/A would 
necessitate. Systems of differential equations which exhibit this behaviour are generally 
referred to as stiff systems. We refrain from formulating a rigorous definition of stiffness. 

4.1 Stability of numerical methods for stiff systems 

In order to motivate the various definitions of stability which occur in this section, we begin 
with a simple example. Consider Euler's implicit method for the initial value problem 

y' = ^y , y(o) = yo , 
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Figure 2: The functions u and v plotted against x for x G [0, 1]. 

where A is a complex number. The stability polynomial of the method is tt{z, h) = p{z) — 
ha{z) where h = hX, p{z) = z — 1 and a{z) = z. Since the only root of the stability 
polynomial is z = 1/(1 — /i), we deduce that the method has the region of stability 

n = {h : \l-h\> 1} . 

In particular TZ includes the whole of the left-hand complex half plane. The next definition 
is due to Dahlquist (1963). 

Definition 11 A linear multistep method is said to be A-stable if its region of absolute 
stability, TZa contains the whole of the left-hand complex half-plane 5R(/iA) < 0. 

As the next theorem by Dahlquist (1963) shows, Definition 11 is far too restrictive. 

Theorem 16 

(i) No explicit linear multistep method is A-stable. 

(a) The order of an A-stable implicit linear multistep method cannot exceed 2. 

(Hi) The second-order A-stable linear multistep method with smallest error constant is 
the trapezium rule. 

Thus we adopt the following definition due to Widlund (1967). 
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Definition 12 A linear multistep method is said to be A{a)-stable, a G (0, 7r/2), if its 
region of absolute stability TZa contains the infinite wedge 

Wa = {h\TT — a < arg(^) < it + a} . 

A linear multistep method is said to be A(0) -stable if it is A{a) -stable for some a G 
(0,7r/2). A linear multistep method is Aq stable if TZa includes the negative real axis in 
the complex plane. 

Let us note in connection with this definition that if 3ftA < for a given A then 
h = hX either Hes inside the wedge Wa or outside Wa for all positive h. Consequently, 
if all eigenvalues A of the matrix A happen to lie in some wedge Wp then an A(/3)-stable 
method can be used for the numerical solution of the initial value problem (85) without 
any restrictions on the step size h. In particular, if all eigenvalues of A are real and 
nonnegative, then an ^(0) stable method can be used. The next theorem (stated here 
without proof) can be regarded the analogue of Theorem 16 for the case of A{a) and ^(0) 
stability. 

Theorem 17 

(i) No explicit linear multistep method is A(0) -stable. 

(ii) The only A{0) -stable linear k-step method whose order exceeds k is the trapezium 
rule. 

(Hi) For each a £ [0,7r/2) there exist A{a)-stable linear k-step methods of order p for 
which k = p = S and k = p = A. 

A different way of loosening the concept of ^-stability was proposed by Gear (1969). 
The motivation behind it is the fact that for a typical stiff problem the eigenvalues of the 
matrix A which produce the fast transients all lie to the left of a line = —a, a > 0, in 
the complex plane, while those that are responsible for the slow transients are clustered 
around zero. 

Definition 13 A linear multistep method is said to be stiffly stable if there exist positive 
real numbers a and c such that TZa ^ TZi U TZ2 where 

TZi = {heC : ^h< -a} and TZ2 = {h e C : -a < m < 0, -c<Qh<c} . 

It is clear that stiff stability implies yl(a)-stability with a = arctan(c/a). More gener- 
ally, we have the following chain of implications: 

A-stability ^ stiff-stability ^ yl(a)-stability =^ ^(O)-stability ^ ^o-stability . 

In the next two sections we shall consider linear multistep methods which are par- 
ticularly well suited for the numerical solution of stiff systems of ordinary differential 
equations. 
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Table 3: Coefficients, order, error constant and stability parameters for backward differ- 
entiation methods 



4.2 Backward differentiation methods for stiff systems 

Consider a linear multistep method with stability polynomial 'K{z,h) = p{z) — ha{z). If 
the method is ^(a)-stable or stiffly stable, the roots r(h) of 7r(-, h) lie in the closed unit 
disk when h is real and h ^ —00. Hence, 

= _lim ^i!lM = _iini a{r{h)) = aiJira r{h)) ; 

h—*—oo h ft— »— 00 ft.— »— 00 

in other words, the roots of 7r(-, h) approach those of a{-). Thus it is natural to choose a{-) 
in such a way that its roots lie within the unit disk. Indeed, a particularly simple choice 
would be to take a{z) = I3kz^] the resulting class of, so-called, backward differentiation 
methods has the general form: 

n 

where the coefficients aj and Pk are given in Table 3 which also displays the value of a in 
the definition of stiff stability and the angle a from the definition of A(a) stability, the 
order p of the method and the corresponding error constant Cp+i for p = 1, . . . , 6. For 
p > 7 backward differentiation methods of order p of the kind considered here are not 
zero-stable and are therefore irrelevant from the practical point of view. 



57 



4.3 Gear's method 



Since backward differentiation metfiods are implicit, tfiey fiave to be used in conjunction 
with a predictor. Instead of iterating the corrector to convergence via a fixed point itera- 
tion, Newton's method can be used to accelerate the iterative convergence of the corrector. 
Rewriting the resulting predictor-corrector multi-step pair as a one step method gives rise 
to Gear's method which allows the local alteration of the order of the method as well 
as of the mesh size. We elaborate on this below. 

As we have seen in Section 4.1, in the numerical solution of stiff systems of ordinary 
differential equations, the stability considerations highlighted in parts (i) of Theorems 16 
and 17 necessitate the use of implicit methods. Indeed, if a predictor-corrector method is 
used with a backward differentiation formula as corrector, a system of nonlinear equations 
of the form 

Yn+k - hl3kiiXn+k,yn+k) = Sn+k 

will have to be solved at each step, where 

fe-i 

3=0 

is a term that involves information which has already been computed at previous steps 
and can be considered known. If this equation is solved by a fixed-point iteration, the 
Contraction Mapping Theorem will require that 

Lh\h\ < 1 (89) 

in order to ensure convergence of the iteration; here L is the Lipschitz constant of the 
function f (x, •). In fact, since the function f (x, •) is assumed to be continuously differen- 
tiable. 



L = max 

(rr,y)eR 



For a stiff system L is typically very large, thus the restriction on the steplength h expressed 
by (89) is just as severe as the condition on h that one encounters when using an explicit 
method with a bounded region of absolute stability. In order to overcome this difficulty. 
Gear proposed to use Newton's method: 



^n+k y n+k 



I - h/3kQ^{Xn+k,yn+k) [yn+k-hPki{Xn+k,yn+k)-Sn+k] , (90) 



for s = 0, 1, . . . , with a suitable initial guess yl^\_i^- Even when applied to a stiff system, 
convergence of the Newton iteration (90) can be attained without further restrictions on 
the mesh size h provided that we can supply a sufficiently accurate initial guess y^^+k (^^ 
using an appropriately accurate predictor, for example). 

On the other hand, the use of Newton's method in this context has the disadvantage 
that the Jacobi matrix df/dy has to be reevaluated and the matrix / — /i/?fc^(x„+fc; Vn+k) 
inverted at each step of the iteration and at each mesh point Xn+k- 

One aspect of Gear's method is that the matrix I — h(3k^{xn+ki Vn+k) involved in the 
Newton iteration described above is only calculated occasionally (i.e. at the start of the 
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iteration, for s = 0, and thereafter only if the Newton iteration exhibits slow convergence) ; 
the inversion of this matrix is performed by an LU decomposition. 

A further aspect of Gear's method is a strategy for varying the order of the backward 
differentiation formula and the step size according to the intermediate results in the com- 
putation. This is achieved by rewriting the multistep predictor-corrector pair as a one-step 
method (in the so-called Nordsieck form). For further details, we refer to Chapter III. 6 in 
the book of Hairer, Norsett and Wanner. 



5 Nonlinear Stability 

All notions of stability which were considered so far in these notes rest on the assumption 
that f{x,y) = Xy or f(a;,y) = Ay + b. The purpose of this section is to develop an 
appropriate theoretical framework which is directly applicable to the stability analysis of 
numerical methods for nonlinear ODEs. 

Consider the linear system of differential equations y' = Ay, y(xo) = yo, where all 
eigenvalues of the m x m matrix A have negative real part. Then ||y(a;)|| decreases as x 
increases; also, neighbouring solution curves get closer as x increases: if u(x) and v(x) are 
two solutions to y' = subject to u(xo) = Uq and v(xo) = Vq, respectively, then 

||u(x2) - v(x2)|| < ei^^-^i)^^i^^M)\\u{xi) - v(a;i)|| , X2 > xi > xq , (91) 

where || • || denotes the Euclidean norm on R"*. Clearly, 

< g(=^2-a;i)maxj 3ftAi(A) ^ ^ 

for X2 > xi > Xq and therefore, 

||U(X2) - V(X2)|| < ||u(xi) - v(xi)|| , X2 > Xl > Xo . 

It is the property (91) that has a natural extension to nonlinear differential equations and 
leads to the following definition. 

Definition 14 Suppose that u and v are two solutions of the differential equation y' = 
f(x,y) subject to respective initial conditions u(xo) = Uq, v(xo) = Vq. // 

||u(x2) - v(a;i)|| < ||u(a;i) - v(a;i)|| 

for all real numbers x^ and xi such that X2 > xi > xq, where \\ ■ \\ denotes the Euclidean 
norm on C^, then the solutions u and v are said to be contractive in the norm || • ||. 

To see what assumptions of f we must make to ensure that two solutions u and v 
to the differential equations y' = f{x,y), with respective initial conditions u{xo) = Uq, 
v(xo) = Vq, are contractive, let (•, ) denote the inner product of R"* and consider the 
inner product of u' — v' = i{x, u) — i{x, v) with u — v. This yields, 

^^||u(x) - v(x)f = (f(x,u(a;)) - f (x, v(x)), u(x) - v(a;)) . 

Thus, if 

(f (x, u(x)) - i{x, v(x)), u(x) - v(x)) < (92) 
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for all X > xq then 

^ '|u(x)-v(x)f <0 



2dx' 

for all X > xq, and therefore the solutions u and v corresponding to initial conditions uo 
and vq, respectively, are contractive in the norm || • ||. The inequality (92) motivates the 
following definition. 

Definition 15 The system of differential equations y' = f{x,y) is said to be dissipative 

in the interval [xo,oo) if 

(f(a;,u) -f(a;,v),u- v) < (93) 
for all X > xo and all u and v in R"^ . 

Thus we have proved that if the system of differential equations is dissipative then any 
solutions u and v corresponding to respective initial values Uq and Vq are contractive. A 
slight generalisation of (92) results in the following definition. 

Definition 16 The function f{x,-) is said to satisfy a one-sided Lipschitz condition 

on the interval [a;o,oo) if there exists a function i/(x) such that 

(f(x,u) -f(a;,v),u- v) < z/(x)||u - (94) 

for all X G [xq, oo). 

In particular, if f (x, •) satisfies a one-sided Lipschitz condition on [xq, oo) and z/(x) < for 
all x G [a;o,oo), then the differential equation y' = f(x,y) is dissipative on this interval, 
and therefore any pair of solutions u and v to this equations, corresponding to respective 
initial values Uq and vq are also contractive. 

Now we shall consider numerical methods for the solution of an initial value problem 
for a dissipative differential equation y' = f(a;,y) and develop conditions under which nu- 
merical solutions are also contractive in a suitable norm. In order to keep the presentation 
simple we shall suppose that f is independent of x and, instead of a linear A;-step method, 

k k 

Yl ^jyn+j =hY^ PAYn+j) (95) 
j=0 j=0 

for the numerical solution of y' = f (y), y(xo) = yo, we shall consider its one- leg twin 

k f \ 

^jyn+j = hfiJ2 (^jyn+j ■ (96) 

i=0 \j=0 J 

For example, the one-leg twin of the trapezium rule method 

Yn+i -yn = ^h [f (y„+i) + f (yn)] 
is the implicit midpoint rule method 

yn+i - yn = /if Q(yn+i + yn. 
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Let us recall the notation from Section 3.1 to simplify writing: putting yn+i = Ey^, we 
can write the linear A;-step method (95) as 

p{E)yn = ha{E){{yn) , 

where p{z) = J2j=o '^j^'^ '^i^) = Sj=o Pj^'' first and second characteristic 

polynomial of the method; the associated one-leg twin (96) is then 

piE)yn = Mia{E)yn) . 

There is a close relationship between the linear multistep method (95) and its one-leg 
twin (96). Let z„ = cj(i?)y„; then 

piE)zn = p{E)a{E)yn = a{E)p{E)yn = ha{E)i{a{E)yn) = /ia(£;)f (z„) . 

In other words, if {yn}n>o is a solution to (96) then {z„}„>o, with z„ = (j(£')y„, is the 
solution of the linear multistep method (95) whose one-leg twin (96) is. This connection 
allows results for the one-leg twin to be translated into results for the linear multistep 
method. 

Now we shall state a definition of nonlinear stability due to Dahlquist (1975). Before 
we do so, we introduce the concept of G-norm. Consider a vector 

— (,z„^j._i, . . . ,z„; 

in R"*'^ where Zn+j G R™, j = 0, 1, . . . , A; — 1. Given that G = (gij) is & k x k sjmimetric 
positive definite matrix, the G-norm || • \\g is defined by 

/ k k \ 

llZnllc = I 9ij{^n+k-i: ^n+k-j) 

\i=lj=l 

Definition 17 The k-step method (96) is said to he G-stable if there exists a symmetric 
positive definite matrix G such that 

||Z„+i|||- ||Z„||| < 2{p{E)zn,a{E)zr,)/a''{\) . 

Let {u„} and {v^} be two solutions of the differential equation y' = f(y) given by 
(96) with different starting values, and suppose that (96) is G-stable. Define the vectors 
U„ and Vn in R"*'^ by 

Un = (u^+A;_i, . . . ,U^)'^ , Vn = (v^^^j.^, . . . , V^)"^ . 

Since we are assuming that the differential equation y' = f(y) is dissipative, it follows 
that 

||Un+i - Vn+i||^ - ||Un - YnWl < '^{p{E){un - v^), a(£;)(u„ - Vn))/a2(l) 

< 2{M{a{E)nn) - /if (a(£;)v„), a(E)(u„ - Vn))/a2(l) 
2h 

= ^^(f (a(^)un) - M{a{E)wn),a{E)Mn - a{E)wn)) < . 
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Hence, 

||U„+i - V„+i||^ < ||U„ - V„||^ forn > . 

Thus wc deduce that a G-stable approximation of a dissipative ordinary differential equa- 
tion y' = f (y) is contractive in the G-norm; this is a discrete counterpart of the property 
we estabUshed at the beginning of this section that analytical solutions to a dissipative 
equation of this kind are contractive in the Euclidean norm. For further developments of 
these ideas, we refer to the books of J.D. Lambert, and A.M. Stuart and A.R. Humphries. 

6 Boundary value problems 

In many applications a system of m simultaneous first-order ordinary differential equations 
in m unknowns yi{x) , y2{x) , . . . , ymix) has to be solved. If each of these variables satisfies 
a given condition at the same value a of 2; then we have an initial value problem for a 
system of first-order ordinary differential equations. If the yi, i = 1, . . . , m, satisfy given 
conditions at different values a, b, c, . . . of the independent variable x then we have a 
multi-point boundary value problem. In particular, if conditions on the yi, i = 1, . . . , m, 
are imposed at two different values a and b then we have a two-point boundary value 
problem. 

Example 7 Here is an example of a multipoint (in this case, three-point) boundary value 
problem: 

y"'-y" + y'-y = 0, y(0) = 1 , y{l) = e, y'{2) = e^. 
The exact solution is y{x) = e^. 

Example 8 This is an example of a two-point boundary value problem: 

y"-2y' = 0, y(l) = l, y'{2) + [y{2)f = 0. 

The exact solution is y{x) = l/x. 

In this section we shall consider three classes of methods for the numerical solution of 
two-point boundary value problems: shooting methods, matrix methods and collocation 
methods. 

6.1 Shooting methods 

Let us consider the two-point boundary value problem 

y" = f{x,y,y') , y{a) = A , y{b) = B , (97) 

with a <b and x G [a, b]. We shall suppose that (97) has a unique solution. The motivation 
behind shooting methods is to convert the two-point boundary value problem into solving 
a sequence of initial value problems whose solutions converge to that of the boundary 
value problem, so that one can use existing software developed for the numerical solution 
of initial value problems: observe that an attempt to solve the boundary value problem 
(97) directly will lead to a coupled system of nonlinear equations whose solution may be 
a hard problem. 
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Let us make an initial guess s for y'{a) and denote by y{x; s) the solution of the initial 
value problem 

y" = f{x,y,y') , y{a) = A , y'{a) = s . (98) 

Introducing the notation u{x;s) = y{x;s), v{x;s) = ^y{x;s), we can rewrite (98) as a 
system of first-order ordinary differential equations: 

d 

7— ^(x: s) = v(x:s), u(a:s)=A, 
ox 

(99) 

d 

— w(x;s) = f{x,u{x;s),v{x;s)), v{a;s) = s. 

The solution u{x; s) of the initial value problem (99) will coincide with with the solution 
y{x) of the boundary value problem (97) provided that that we can find a value of s such 
that 

(f){s) = u{b;s)-B = 0. (100) 

The essence of the shooting method for the numerical solution of the boundary value 
problem (97) is to find a root to the equation (100). Any standard root-finding technique 
can be used; here we shall consider two: bisection of the interval which is known to contain 
the root and Newton's method. 

6.1.1 The method of bisection 

Let us suppose that that two numbers si and S2 are known such that 

(f){si) < and ^(§2) > . 

We assume, for the sake of definiteness, that si < S2- Given that the solution of the initial 
value problem (99) depends continuously on the initial data, there must exist at least one 
value of s in the interval (si,S2) such that ^(s) = 0. Thus the interval [si,S2] contains a 
root of the equation (100). 

The root of (100) can now be calculated approximately using the method of bisection. 
We take the midpoint S3 of the interval [si,S2], compute u(6, S3) and consider whether 
0(53) = u{b] S3) — B is positive or negative. If ^(53) > then it is the interval [51,53] 
that contains a root of 0, whereas if ^(53) < then the interval in question is [53, 52]- By 
repeating this process, one can construct a sequence of numbers {sn}^=i converging the 
s. In practice the process of bisection is terminated after a finite number of steps when 
the length of the interval containing s has become sufficiently small. 

6.1.2 The Newton— Raphson method 

An alternative to the method of bisection is to compute a sequence {sn}^=i generated by 
the Newton-Raphson method: 

Sn+l = Sri - 0(Sn)/0'(-Sn) , (101) 

with the starting value sq chosen arbitrarily in a sufficiently small interval surrounding the 
root. For example, a suitable sq may be found by performing a few steps of the method of 
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bisection. If sq is a sufficiently good approximation to tfie required root of (100) the theory 
of the Newton-Raphson method ensures that, in general, we have quadratic convergence 
of the sequence {sn}^=Q to the root s. 

Prom the point of view of implementing (101) the first question that we need to clarify 
is how one can compute compute 0'(s„). To do so, we introduce the new dependent 
variables 

C(.;.) = ^, .(x;.) = ^ 

and differentiate the initial value problem (99) with respect to s to obtain a second initial 
value problem: 



d^jx; s) 
dx 

dri{x; s) 
dx 



where 



ri{x;s), ^(o;s) = 0, 

p{x; s)^{x; s) + q{x; s)r]{x; s) , r/(a; s) = 1 , 
df{x,u{x;s),v{x;s)) 



(102) 



p{x;s) 
q{x;s) 



du 

df{x,u{x;s),v{x;s)) 



(103) 



dv 

If we assign the value s„ to s, n > 0, then the initial value problem (99), (102) can be solved 
by a predictor-corrector method or a Runge-Kutta method on the interval [a, b]. Thus we 
obtain u{b] Sn) or, more precisely, an approximation to u{b; Sn) from which we can calculate 
= u{b; Sn) — B; in addition, we obtain an approximation to ^{b; s„) = (f)'{sn)- Having 
calculated (pisn) and 4>'{sn), we obtain the next Newton-Raphson iterate from (101). 
The process is then repeated until the iterates Sn settle to a fixed number of digits. 
Two remarks arc in order: 

1) According to (103), the initial value problems (99) and (102) are coupled and there- 
fore they must be solved simultaneously over the interval [a,b], with s set to Sn, 
n = 0,l,2,...; 

2) The coupled initial value problem (99), (102) may be very sensitive to perturbations 
of the initial guess so; a bad initial guess of so may result in a sequence of Newton- 
Raphson iterates {s„}^q which does not converge to the root s. 

The latter difficulty may be overcome by using the multiple shooting method 
which we describe below. First, however, we show how the simple shooting method may 
be extended to the nonlinear two-point boundary value problem 

y' = f(x,y), a<x<b, g(y(a), y(6)) = , (104) 

where y, f and g are m-component vector functions of their respective arguments. In the 
simple shooting method the boundary value problem (104) is transformed into the initial 
value problem 

u'(x; s) = f(x, u(x; s)) , a<x<b, u(o;s) = s, (105) 
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whose solution is required to satisfy the condition 

G(s) = g(s,u(6;s)) = (106) 

for an unknown value of s. Thus the problem has been transformed into one of finding a 
solution to the system of nonlinear equations G(s) = 0. In order to evaluate the function 
G for a specific value of s a numerical method for the solution of the initial value problem 
(105) has to be employed on the interval [a,b] to compute (an approximation to) u(6;s). 
In fact, as we have already seen earlier, a root-finding procedure such as Newton's method 
will require the computation of the Jacobi matrix 

This in turn requires the solution of a coupled initial value problem for linear ordinary 
differential equations to obtain du/ds for a < x < b. 

The shooting method is said to converge if the root-finding algorithm results in a 
sequence {s„}^q which converges to a root s of G; then s = y(a) where y(a;) is the 
desired solution of the boundary value problem (104). The convergence of this sequence 
and therefore the success of the shooting method may be hampered by two effects: 

1) A well-conditioned boundary value problem of the form (104) may easily lead to an 
ill-posed initial value problem (104); 

2) A bounded solution to the initial value problem (104) may exist only for s in a small 
neighbourhood of y(s) (which, of course, in unknown). 

The idea behind the multiple shooting method is to divide the interval [a, b] into 
smaller subintervals 

a = xo < xi < . . . < xk-1 < xk = b ; (107) 

the problem then is to find a vector S"^ = (sq , • • • , sj-.^) such that the solutions Uk{x; Sk) 
of the initial value problems 

Ufc(^;sfc) = f(x,Ufc(a;;Sfc)) , Xk < x < Xk+i , 

(108) 

Ufc(xjk;Sjk) = Sjfc , k = 0,...,K -1 , 

satisfy the conditions 

Ufe(a;jk+i;Sjt) - Sjfc+i = , k = 0,...,K -2 , g{so,UK-i{b;SK-i)) = ■ (109) 

The equations (109) can be written in the compact form G(S) = 0. A clear advantage 
of the multiple shooting method over simple shooting is that the growth of the solutions 
to the initial value problems (108) and the related linear initial value problems for the 
duk{x; Sk)/dsh, k = 0,. . ., K—1, can be approximated accurately by selecting a sufficiently 
fine subdivision (107) of the interval [a,b]. 

Indeed, it is possible to prove that, under reasonable assumptions, the multiple shoot- 
ing method leads to an exponential increase of the size of the domain of initial values 
for which the first iteration of the root-finding procedure is defined. This is consistent 
with the practical observation that multiple shooting is less sensitive to the choice of the 
starting values than simple shooting. 
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6.2 Matrix methods 



In this section, rather than attempting to convert the boundary value problem to an 
initial value problem, we approximate it directly by using a finite difference method. This 
results in a system of equations for the unknown values of the numerical solution at the 
mesh points. We begin by considering a linear boundary value problem. In this case 
the calculation of the numerical solution amounts to solving a system of linear equations 
with a sparse matrix. We then consider a nonlinear boundary value problem; upon the 
application of Newton's method, this again involves, in each Newton iteration, the solution 
of a system of linear equations with a sparse matrix. Methods of this kind are usually 
referred to in the ODE literature as matrix methods. 



6.2.1 Linear boundary value problem 

Let us consider the two-point boundary value problem 

y" +p{x)y' + q{x)y = f{x) , x e {a,b) , 

aoy{a) + boy'{a) = cq , (110) 
aiy{b) + biy'{b) = ci . 

We discretise (110) using a finite difference method on the uniform mesh 

{xi \ Xi = a + ih , i = 0, . . . , N} 

of step size h = {b — a)/N, N > 2. The essence of the method is to approximate the 
derivatives in the differential equation and the boundary conditions by divided differences. 

Assuming that y is sufficiently smooth, it is a simple matter to show using Taylor series 
expansions that 

y"ix,) = y(-^^^)-^y(f+y(-^-^) +o{h^), (m) 

y'ix.) = ^(^m)-j(x.-i) ^^(^.^^ ^^^2) 

y'ixo) = zM^^^l±^M^yM + o{h'), (113) 
y'{x^) = y±I^^^l^Ml^^ll±M^ + o{h^). (114) 

Now, using (111-114) we can construct a finite difference method for the numerical solution 
of (110); we denote by j/j the numerical approximation to y{xi) for i = 0, . . . , N: 

Vi+i '^Vi + Vi-i = f{xi), i = l,...,N -1 , 



-3yo + 4yi - y2 

aoyo + oo ^1 = Co, 

in 

, , yN-2 - i^A -i + 3yN 

aiyN + bi — = ci . 

2h 

After rearranging these, we obtain 



66 



bi 46i / 35i 

This is a system of linear equations of the form 

^oVo + Coyi + Boy2 = Fq , 
Aiyi^i + dyi + Biyi+i = Fi , 
ANyN-2 + CNyN-i + B^yN = Fn 

The matrix of the system is 



bo 

yN-2 



Co , 

Ci . 



l,...,iV-l , 



r ^0 


Co 


Bo 












Bi 











A2 


C2 


B2 

















An-1 














An 






. . ... 

Cn-1 Bn-i 
Cn Bn 

Let us consider the following cases: 

1) If So = and ^jv = then M is a tridiagonal matrix. 

2) If Bo 7^ and Bi = (and/or A^ / and ^at-i = 0) then we can interchange 
the first two rows of the matrix (and/or the last two rows) and, again, we obtain a 
tridiagonal matrix. 

3) If -Bq 7^ and Bi ^ (and/or A^ ^ and ^at-i 7^ 0) then wc can eliminate Bq 
from the first row (and/or An from the last row) to obtain a tridiagonal matrix. 

In any case the matrix M can be transformed into a tridiagonal matrix. Thus, from now 
on, without any restrictions on generality, we shall suppose that M is tridiagonal. To 
summarise the situation, we wish to solve the system of linear equations 



My = F 



where M is a tridiaginal matrix, 

y = {yo,yi,---,yN)'^, 



F = {Fo,Fi, . . . ,Fn)'^ 



The algorithm we present below for the solution of this linear system is usually referred 
to as the Thomas algorithm. It is a special case of Gaussian elimination or LU 
decomposition; in fact, it is this latter interpretation that we adopt here. 
We wish to express M as a product LU of & lower-triangular matrix 



L = 



1 

h 









1 

12 









1 























1 

In 









1 
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and an upper-triangular matrix 



U = 



uo vq 
m 

U2 "^2 









UN-1 VN-1 

UN 



On multiplying L and U and equating the resulting matrix with M, we find that 

i = l,...,N. 



hUi—l — Ai 

uo = Co , vo = Bq , kvi-i +Ui =Ci 

Vi = Bi 



Hence 

Vi = Bi , i = 0,...,N , no = Co , 



Ui =Ci- {AiBi^i)/ui-i 
h — Ai/ui—i 



i = l,...,N. 



Given that the entries of the matrix M are known, the components of L and U can be 
computed from these formulae, and the set of linear equations My = F can be written in 
the following equivalent form: 

f Lz = F 
[Uy = z ■ 

Thus, instead of solving My = F directly, we solve two linear systems in succession, each 
with a triangular matrix: Lz = F is solved for z, followed by solving ?7y = z for y. 
Writing this out in detail yields 



and 



li—lZi 



VN = zn/un , 
Vi = {zi - Viyi+i)/ui 



i = N -1,...,0 . 

Expressing in these formulae the values of Ui and Vi in terms of Ai, Bi, Ci, we find that 

,0, 



Vi = ai+iyi+i + (3i+i , i = N -1, 
VN = Pn+1 , 



where 



Bi 



Ci + OiAi 
Fi - PiAi 



1,2, 



1,2,. ..,N 



ai 

Pi 



Co ' 
Fo 
Co ■ 



Ci + aiAi ' 

The last set of formulae are usually referred to as the Thomas algorithm. 
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6.2.2 Nonlinear boundary value problem 

Now consider a nonlinear second-order differential equation subject to linear boundary 
conditions: 

y" = f{x,y,y'), xe{a,b), 

aoy{a) + boy'{a) = cq , 
aiy{b) + biy'{b) = a . 

On approximating the derivatives with divided differences, we obtain the following set of 
difference equations: 

ym-2..+..-i _ f(.^,yJ-^±^^] , ^=l,...,N-l, 



-3yo + Ayi - y2 

aoyo + oo — = Co, 

2h 

aiyN + h — = ci . 

2h 

After rearranging these, we obtain 

J^Vi-i ' J^Vi + J^Vi+i = f(xi,yi, ^ ), l<i<N-l 

36o\ , 46o bo 



bi Abi ( 36i 



This is a system of system of nonlinear equations of the form 

Aiyi-x + Ciyi + B^yi^x = Fi , i = l,...,N -1 , (115) 

ANyN-2 + CNyN-l + BjsfyN = , 

where 

Ao = ao- (36o)/(2/i) , Co = 4bo/{2h) , Bo = -bo/{2h) , 

Ai = l/h? , d = -2/h\ Bi = l//i2, z = l,...,iV-l , 

An = bi/{2h), Cn = -46i/(2/i), B^ = ai + (36i)/(2/t) , 

and 

Fo = Co , Fi = f{xi, yi, {yi+i - yi-i)/ {2h)) , i = l,...,N -I , Fn = ci . 
Thus, (115) can be written in the compact form 

My = F(y) 

where M is a tridiagonal matrix, y = (yo, • • • , Z/at)^ and F(y) = (Fo, . . . , F^)^ ■ The 
nonlinear system of equations 

G(y) = My - F(y) = 
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can now be solved by Newton's method, for example, assuming that a solution exists. 
Given a starting value y'^'^^ for the Newton iteration, subsequent iterates are computed 
successively from 

j(yH)(y(n+i) _ y(n)) ^ _G(yH) , n = 0, 1, 2, . . . , (116) 

where 

\^yj Jo<i,j<N 

is the Jacobi matrix of G. As Gi is independent of yj with \i — j\ > 1, the matrix J(y("^) 
is tridiagonal. Consequently, each step of the Newton iteration (116) involves the solution 
of a system of linear equations with a tridiagonal matrix; this may be accomplished by 
using the Thomas algorithm described above. 



6.3 Collocation method 

Consider the boundary value problem 

= y" + Pix)y' + q{x)y = f{x) , x G (a, 6) , 

aoy{a) + hoy' {a) = cq , (117) 
aiy{h) + hiy'{h) = ci , 

where oq, ai, 6o, hi are real numbers such that 

(ao&i — ai6o) + aoai(6 — a) 7^ . (US) 

It can be assumed, without loss of generality, that cq = and ci = 0; for if this is not the 
case then we consider the function y defined by 

y{x) = y{x) - d{x) , 

with d{x) = ax + (3 where a and (3 are real numbers chosen so that d{x) satisfies the 
(nonhomogeneous) boundary conditions at .x = a and x = b stated in (117). It is a 
straightforward matter to see that provided (118) holds there are unique such a and p. 
The function y then obeys the differential equation Cy = f, where f{x) = f{x) — Cd{x), 
and satisfies the boundary conditions in (117) with cq = ci = 0. 

Suppose that {ipi}fLi is a set of linearly independent functions defined on the interval 
[a,h] satisfying the homogeneous counterparts of the boundary conditions in (117) (i.e. 
Co = ci = 0). We shall suppose that each function i/ji is twice continuously differentiable 
on [a, h]. 

The essence of the collocation method is to seek an approximate solution yN{x) to 
(110) in the form 

N 

Vn{x) = Y,iii^i{x) , (119) 

i=l 

and demand that 

CyN{xj) = f{xj), j = l,...,N, (120) 
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at (N+l) distinct points xj, j = I, . . . , N, referred to as the collocation points. We note 
that since each of the functions ipi{x) satisfies the (homogeneous) boundary conditions at 
X = a and x = b the same is true of yN{x)- 

Now, (119-120) yield the system of linear equations 

N 

^ijCAixj) = fixj) , j = l,...,N, (121) 

1=1 

for the coefficients ^i, i = 1,...,N. The specific properties of the collocation method 
depend on the choice of the basis functions ipi and the collocation points xj. In general 
the matrix M = {C^i{xj))i<ij<j\f is full. However, if each of the basis functions Vi has 
compact support contained in (a, b) (for example, they are B-splines) then M is a band- 
matrix, given that both it'iixj) = and Cipi{xj) = for |i — j| > M for some integer M, 
1 < M < A^. If M is a fixed integer, independent of iV, then the resulting system of linear 
equations can be solved in 0{N) arithmetic operations. 

In spectral collocation methods, the functions il^i{x) are chosen as trigonometric poly- 
nomials. Consider, for example, the simple boundary value problem 

-y"ix) + qix)y{x) = fix) , xG(0,7r), y(0) = y(7r) = , 

where q{x) > for all x in [0, vr]. The functions 'ipiix) = sin ix, i = 1, . . . ,N , satisfy the 
boundary conditions and are linearly independent on [0, tt]. Thus we seek an approximate 
solution yN{x) in the form 

N 

yN{x) = Y^isinix . 

i=l 

Substitution of this expansion into the differential equation results in the system of linear 
equations 

N 

Yl 0^ + (ii^j)) sinixj- = f{xj) , i = l,...,N, 

i=l 

for ^1, . . . ,^jv. The collocation points are usually chosen as 

Xo = , 7 = 1, . . . , TV . 

This choice is particularly convenient when g is a (nonnegative) constant, for then the 
linear system 

N 

^ei(^' + g)sin-^ = /(x,) , i = l,...,N, 

can be solved for ^i, . . . , in 0(N) arithmetic operations using a Fast Fourier Transform, 
despite the fact that the matrix of the system is full. 
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Further Exercises 



1. Verify that the following functions satisfy a Lipschitz condition on the respective 
intervals and find the associated Lipschitz constants: 



a) f{x,y) = 2yx'^ , x G [l,oo) ; 

b) fix, y) = e-=^^ tan-^ y , x G [1, oo) ; 

c) f{x, y) = 2y{l + y'^)-^{l + e"^) , x G (-00, cx)) . 

2. Suppose that m is a fixed positive integer. Show that the initial value problem 

y' = y2W(2m+l) ^ y(o) ^ ^ 

has infinitely many continuously differentiable solutions. Why does this not contra- 
dict Picard's theorem? 

3. Show that the explicit Euler method fails to approximate the solution y{x) = 
(4x/5)^^^ of the initial value problem y' = y^^^, y{0) = 0. Justify your answer. 
Consider the same problem with the implicit Euler method. 

4. Write down the explicit Euler method for the numerical solution of the initial value 
problem y' + 5y = xe~^^ , y(0) = , on the interval [0, 1] with step size h = 1/N, 
N > 1. Denoting by y^ the Euler approximation to y{l) at x = 1, show that 
lim^v^ooyiv = Find an integer Nq such that 

|y(l)-yiv| < 10-^ foraUAr>iVo. 

5. Consider the initial value problem 

y' = log log(4 + y^), x G [0, 1] , y(0) = 1 , 

and the sequence {yn}n=o^ N > 1, generated by the explicit Euler method 

yn+i = yn + ^loglog(4-|-y^) , n = 0,...,N -1 , yo = 1 , 

using the mesh points x„ = nh, n = 0,...,N, with spacing h = 1/N. Here log 
denotes the logarithm with base e. 

a) Let Tn denote the truncation error of Euler's method at x = x„ for this initial 
value problem. Show that |T„| < h/4. 

b) Verify that 

|y(x„+i) - y„+i| < (1 + hL)\y{xn) - yn\ + h\Tn\ , n = 0, ...,A^-1, 

where L = 1/(2 log 4). 

c) Find a positive integer Nq, as small as possible, such that 



max \v(xn) — Vn\ < 10 



whenever N > Nq. 
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6. Define the truncation error T„ of the trapezium rule method 



Vn+l =yn + ^h {fn+1 + fn) 

for the numerical solution of the initial value problem y' = f{x, y), y(0) given, where 

fn = f{Xn, Vn) and h = Xn+l - Xn- 

By integrating by parts the integral 



/ {x - Xn+l){x - Xn)y"'ix)dx , 

J Xn. 



or otherwise, show that 



1 

12 



for some in the interval (xjj, Xn+i), where y is the solution of the initial value 
problem. 

Suppose that / satisfies the Lipschitz condition 

\f{x,u) - f{x,v)\ < L\u-v\ 

for all real x, u, v, where L is a positive constant independent of .x, and that |y'"(a;)| < 
M for some positive constant M independent of x. Show that the global error 
Cn = y{xn) — yn Satisfies the inequality 

|en+i| < |en| + ^hL (|e„+i| + \en\) + ^^^^ ■ 
For a uniform step h satisfying hL < 2 deduce that, if yo = y{xo), then 



12L 



1 + \hL" 



7. Consider the following one-step method for the numerical solution of the initial value 
problem y' = f{x, y), y{xo) = yo- 

1, 



Vn+l =yn+ + k2) , 



where 



kl = f{Xn,yn), k2 = f{Xn + h,yn + hki) . 

Show that the method is consistent and has truncation error 



1, 
6' 



fyifx + fyf) c^fxx + 2/12// + jyyj ) 



8. Show that for = 1, 2 there is no i?-stage Runge-Kutta method of order i? + 1. 
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9. Consider the one-step method 

Vn+i = yn + h{aki +bk2) , 

where 

^1 = f{.Xni Vn) ) 

k2 = f{xn + ah, Un + Phki) , 

and where a, b, a, (3 are real parameters. Show that there is a choice of these 
parameters such that the order of the method is 2. Is there a choice of the parameters 
for which the order exceeds 2? 

10. Consider the one-step method 

Vn+l =yn + ahf{Xn, Vn) + Phf{Xn + jh, yn + lhf{Xn, Vn)) , 

where a, (3 and 7 are real parameters. Show that the method is consistent if and 
only \i a + (3 = 1. Show also that the order of the method cannot exceed 2. 

Suppose that a second-order method of the above form is applied to the initial value 
problem y' = —Ay, y(0) = 1, where A is a positive real number. Show that the 
sequence {yn)n>o is bounded if and only if h < |. Show further that, for such A, 

\y{Xn) - Vnl < ^A^/l^Xn , n > . 

11. Derive the mid-point rule method and the Simpson rule method by integrating the 
differential equation y' = f{x,y{x)) over suitable intervals of the real line and ap- 
plying appropriate numerical integration rules. 

12. Write down the general form of a linear multistep method for the numerical solution 
of the initial value problem y' = f{x,y), y{xo) = yo. What does it mean to say that 
such a method is zero-stable! Explain the significance of zero-stability. What is the 
truncation error of such a linear multistep method? 

Determine for what values of the real parameter b the linear multistep method defined 
by the formula 

Vn+Z + (26 - 3)(2/n+2 - 2/n+l) - Vn = hb{fn+2 + fn+l) 

is zero-stable. Show that there exists a value of b for which the order of the method 
is 4. Show further that if this method is zero-stable then its order cannot exceed 2. 

13. Consider the initial value problem y' = f{x,y), y{Q) = yo. Consider attempting to 
solve this problem by the linear multistep method 

ayn-2 + byn-l + yn+l = hf{Xn, Vn) 

on the regular mesh Xn = nh where a and b are constants. 

a) For a certain (unique) choice of a and b, this method is consistent. Find these 
values of a and b and verify that the order of accuracy is 1. 
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b) Although the method is consistent for the choice of a and b from part a) , the 
numerical solution it generates will not, in general, converge to the solution of 
the initial value problem as /i — > 0, because the method is not zero-stable. Show 
that the method is not zero-stable for these a and b, and describe quantitatively 
what the unstable solutions will look like for small h. 

14. Consider the linear two-step method 

yn+2 -yn = ^(/n+2 + 4/„+l -|- fn) ■ 

Show that the method is zero-stable; show further that it is third-order accurate, 
namely, Tn = 0{h^). 

15. Show that the linear three-step method 

llyn+3 + 27yn+2 - 27y„+i - lli/„ = 3/i[/„+3 + 9/^+2 + 9/„+i + /„] 

is sixth order accurate. Find the roots of the first characteristic polynomial and 
deduce that the method is not zero-stable. 

16. Write down the general form of a linear multi-step method for the numerical solution 
of the initial value problem 

y' = f{x, y) , y{xo) = yo , 

on the closed real interval [xq, xj^], where / is a continuous function of its arguments 
and yo is a given real number. Define the truncation error of the method. What 
does it mean to say that the method has order of accuracy p? 

Given that a is a positive real number, consider the linear two-step method 

yn+2 -Oiyn = ^ [f{Xn+2, yn+2) + 4/(Xn+i, -|- , 

on the mesh {xn : Xn = xq + nh,n = 1, . . . , A^} of spacing h, h > 0. Determine 
the set of all a such that the method is zero-stable. Find a such that the order of 
accuracy is as high as possible; is the method convergent for this value of a? 

17. Compute a numerical solution to the initial value problem y' + y = 0, y(0) = 1, on 
the interval [0, 1] with /i = 2"* for A; = 1, 2, . . . , 10, using the linear two-step method 

yn - 2/n-l = ^Hfn^l " fn-2) , 

where the missing starting value yi is computed by the explicit Euler method. Tab- 
ulate the global error at x = 1 for /c = 1, 2, . . . , 10, and comment on its rate of 
decay. 

18. Solve numerically the initial value problem 

y' = x-y\ y{0) = , x G [0, 1] , 

with step size /), = 0.01 by a fourth-order Milne-Simpson method. Use the classical 
fourth-order Runge-Kutta method to compute the necessary starting values. 
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19. Find which of the following linear multistep methods for the solution of the initial 

value problem y' = f(x, y), y{0) given, are zero-stable. For any which are zero-stable, 
find limits on the value of h = Xn+i — Xn for which they are absolutely stable when 
applied to the equation y' = Ay, A < 0. 

a) yn+l -yn = hfn 

b) yn+l +yn- 2yn-l = h{fn+l + fn + fn-l) 
C) yn+l - Vn-l = lHfn+1 + 4/„ -|- fn-l) 

d) yn+l - yn = \h{^fn " /m ) 

e) yn+l -yn = ^h{hfn+l + 8fn - fn-l) 

20. Determine the order of the linear multistep method 

yn+2 - (1 + a)yn+i + yn = [(3 - a)fn+2 + (1 " 3a)/n] 
and investigate its zero-stability and absolute stability. 

21. If (t{z) = is the second characteristic polynomial of a linear multistep method, find 
a quadratic polynomial p(^z) such that the order of the associated linear multistep 
method is 2. Is this method convergent? What is its interval of absolute stability? 

22. Find the interval of absolute stability for the two-step Adams-Bashforth method 

yn+2 - yn+l = ^h [3/n+l - fn] 

using Schur's criterion and the Routh-Hurwitz criterion. 

23. Given that a is a positive real number, consider the linear two-step method 

h 

yn+2 - ayn = [fn+2 + 4/„+l + fn] 

on the mesh {x„ : x„ = + nh, n = 0, . . . . A^} of spacing h, h > 0. For values of 
a such that the method is zero-stable investigate its absolute stability using Schur's 
criterion. 

24. A predictor P and a corrector C are defined by their characteristic polynomials: 

P: p*{z) = z^-l, a*{z) = l{2z^-z^ + 2z) , 
C: p{z) = z^-l, a{z) = liz'^ + 'iz + l) . 

a) Write down algorithms which use P and C in the P{EC)"^E and P{EC)'^ 
modes. 

b) Find the stability polynomials '^p[ec)"^e{z'i h) and TTpi^EC)'" (-^i h) of these meth- 
ods. Assuming that m = 1, use Schur's criterion to calculate the associated 
intervals of absolute stability. 

c) Express the truncation errors Tn^^^^ ^ and Tn^^^^ of these methods in the 
form 0{h^) where r = r{p*,p, m), withp* and p denoting the orders of accuracy 
of P and C, respectively. 
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25. Which of the following would you consider a stiff initial value problem? 

a) y' = -(10^6-^°*^ + l){y-l), y{0) = on the interval x G [0, 1]. Note that the 
solution can be found in closed form: 

y(x) = eio(e-^°'^-i)e- + 1 . 

b) 

y[ = -0.5yi+ 0.5012/2, yi(0) = 1.1 , 
y'^ = 0.5012/1-0.5^2, ^2(0) = -0.9 , 

on the interval x G [0, 10^]. 
c) y' = A{x)y, y(0) = (1, 1, 1)^, where 



A{x) = 



I -1 + 100cos200x 100(1- sin 200x) 

-100(1 + sin 200a;) -(1 + 100cos200x) 

^ 1200(cosl00x + sinl00x) 1200(cos lOOx - sin lOOx) -501 



26. Consider the ^-method 

Vn+l =yn + h[{l- 9)fn + 9fn+l] 

for e e [0,1]. 

a) Show that the method is A-stable for 9 > 1/2. 

b) Show that the method is A(Q;)-stable when it is A-stable. 

27. Show that the second-order backward differentiation method is A-stable. Show that 
the third-order backward differentiation method is not ^-stable, but that it is stiffly 
stable in the sense of Definition 13. 

28. Find Xm > as large as possible such that the system of differential equations 

y'l = -yi + xy2 
y'2 = x^{yi-y2) 

is dissipative in the interval [0,Xm]. Deduce that any two solutions with respective 
initial conditions are then contractive on [0, Xm] in the Euclidean norm. 

29. Show that the trapezium rule method is G-stable with G = 1. 

30. Show that the second-order backward differentiation method and its one-leg twin 
are G-stable with 

^ 5/2 -1 
-1 1/2 



G = 



31. Develop a shooting method for the numerical solution of the boundary value problem 

-/ + ye^ = 1 , y(0) = y(l) = , 
on the interval [0, 1] using: 
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a) The method of bisection; 

b) The Newton-Raphson method. 

Explain how your algorithm can be extended to the case of multiple shooting. 

32. Construct a three-point finite difference scheme for the numerical solution of the 
boundary value problem 

+ x^y = , 2/(0) = 1 , y'{l) = , 

on the interval [0, 1]. Show that the resulting system of equations can be written 
so that its matrix is tridiagonal. Apply the Thomas algorithm to the linear system 
when the spacing between the mesh points \s h = 1/10. 

33. Suppose that real numbers a^, hi and q satisfy 

> , 6j > , Ci>ai + hi, i = l,2,...,N -1 , 

and let 



with eo = 0. Show by induction that < < 1 for i = 1, 2, . . . , — 1, and that the 
conditions 

Cj > , Cj > |aj| + , 1 = 1,2,. .. ,N -1, 
are sufficient for eo = to imply that |ej| < 1 for z = 1, 2, . . . , A/" — 1. 
How is this method used to solve the system of equations 

-aiVi-i + CiUi - biUi+i = fi , i = 1,2, . . . ,N - 1 , 

with 

yo = , yiv = ? 

34. Assume that the boundary value problem 

y" + f{x,y) = 0, 0<x<l, y(0) = y(l) = , 

has a unique solution with a continuous fourth derivative on the interval [0, 1]. Sup- 
pose further that a unique approximate solution can be computed satisfying the 
finite difference scheme 

h~^S^yn + f{Xn, yn) = , l<n<N -1 , yQ = ypf = , 

where 6'^yn = yn+i — 2yn + yn-i, Xn = nh (0 < n < N), and Nh = 1 for some integer 
N>2. 

a) Show that the truncation error of the finite difference scheme is given by 
for some ^„ G (x„_i, x„+i). 
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Show further that the global error e„ = — y„ satisfies 

/i"^(5^e„ + pnCn = Tn , 1 <n < N - 1 , eo = CN = , 

where Pn = fy{xn,f]n) for some r/„ between y{xn) and y„, and it is assumed 
that fy{x, y) is a continuous function of x and y for x € [0, 1] and y G R. 

b) Suppose now that fy{x,y) < for all x £ [0,1] and y G R. Let \Tn\ < M, 
l<n< N-1. Show that Wn = \h?Mn{N - n) satisfies h'^S^Wn = -M, that 

h~'^5'^Wn + PnWn < -M , 1 <n < N - 1 , 

and that, if Vn = Wn + Sn or Vn = Wn — en, then Vn satisfies 

h~^S'^Vn + PnVn <0, 1 < 1% < N — 1 , Vq = Vn = . 

c) Assume that Vn has a negative minimum for some value of n between 1 and 
N — 1 and show that this leads to a contradiction. Deduce that 

vn>0 , l<n<N -1 , 

that 

l^nl ^ Wn , 1 < n < N — 1 , 

and that 

35. Consider the boundary value problem 

-y" + y' + y = x'' , y(0) = , y(l) = . 

Develop a collocation method for the numerical solution of this problem on the 
interval [0, 1] using the collocation points 

7 = 1,..., AT. 
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and the basis functions i/ji{x) = sin(i7rx), i = 1, . . . , N. Solve the resulting system of 
linear equations for increasing values of N and compare the numerical solution with 
the exact solution to the problem. 
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